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ABSTRACT 


Four  techniques  for  the  numerical  solution  of  partial 
differential  equations  and  eigenvalue  problems  were 
investigated.  Typical  problems  considered  were  elliptic 
partial  differential  equations  of  the  form 

u»  +  uyy  =  (1) 

or 

Ux*  *  Uyy  +  (2) 

where  appropriate  boundary  conditions  are  specified  so 
that  the  problem  is  self-adjoint. 

The  four  methods  are  relaxation,  Galerkin,  Rayleigh- 
Ritz,  and  dynamic  programming  combined  with  Stodola’s 
method,  for  eigenvalue  problems. 

The  results  indicated  that  for  eigenvalue  problems 
relaxation  or  dynamic  programming  modified  is  to  be 
preferred  usually  and  for  partial  differential  equations 
Galerkin  or  dynamic  programming  is  preferred. 
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I.  INTRODUCTION 


Initially  for  this  thesis  it  was  planned,  to  investigate 

i  i 

methods  for  finding  eigenfunctions  and  eigenvalues,  with 
a  particular  interest  in  the  oscillation  in  basins  such  as 
harbors  and  bays.  The  report  by  Angel  fio]  introduced 
a  new  mothod,  that  of  dynamic  programming,'  for  the  solu¬ 
tion  of  some  partial  differential  equations.  This  method 
seemed  promising  and  it  was  then  extended  here  to- the  solu¬ 
tion  of  other  partial  differential  equations.  It  also 
made  it  possible  to  invert  a  linear . differential  operator 
and  hence  apply  Stodola's  method  in  finding  eigenvalues 
and  eigenfunctions.  This  led  naturally  to  a  comparison 
with  several  procedures  already  known  to  try  to  compare 

! 

convergence,  speed,  and  accuracy,  by  applying  them  to  the 
solution  of  several  rather  simple  problems. 

It  is  assumed  that  the  reader  has  some  familiarity  with 
the  techniques  of  replacing  differential  equations  by  dif¬ 
ference  equations  and  the  standard  technique  of  separating 
variables  In  linear  partial  differential  equations i  The 
regions  and  their  boundaries  were  assumed  to  be  "nice,'1 
not  excessively  irregular.  Some  knowledge  of  calculus 
of  variations  also  is  desirable  but  not  necessary;  ref¬ 
erence  {5}  provides  more  than  adequate  background. 

Four  problems  were  divided  into  two  categories.  In  the 

i 

first  category  were  eigenvalue  problems.  Two  typical  ones 
were 


6 


u  +-U  *  -  X2U 

xx  yy  A  u 


In  A, 


u.i) 


subject  to  the  constraint 

i 

'  *  0  on  SA,  (1.2) 

where  4A  le  the  boundary  of  the  domain  A;  and  seeond 

?‘U"X!U"  lhA,  (1.3, 

subject  to  the  constraint 

3*U 


u(x.y)  = 


32n2 


*  0 


on  6 A. 


(l.U) 


In  the  second  'category  were  equations  such  as  Poisson’s 
equation  , 


XX  +  Uyy  s  >y ) 


subject  to  the  constraint 
U(x,y)  «  g(x,y ) 
and  the  biharmonic  equation 

I 

s  f(x,y)  : 

subject  to  the  constraints 
U(x,y)  =  h (x ,y ) 

‘  and 


in  A, 


on  6A , 


in  A, 


on  6A 


3  2U  ' 
3n: 


*  q(x,y) 


on  SA. 


(1.5) 


(1.6) 


(1.7) 


(1.8) 


(1.9) 


Problems  in  the  first  category  were  numerically  solved 
by  a  relaxation  method,  the  Rayleigh-Ritz  method,  and 
dynamic  programming  combined  with  Stodola's  method. 
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Problems  in  the  second  category  were  solved  by  Galerkin's 
method,  the  Rayleigh-Ritz  method  and  dynamic  programming. 

The  domain  used  throughout  this  paper  for  comparisons 
was  the  unit  square,  with  the  constraint 

U(x,y)  =0  on  6A .  (1.10) 

Computations  were  also  carried  out  for  L-shaped  and  trian¬ 
gular  regions  and  for  other  boundary  conditions. 

In  the  relaxation  method  and  dynamic  programming  method 
in  solving  eigenvalue  problems  the  initial  estimates  to 
the  eigenfunctions  had  a  special  property.  The  first  ap¬ 
proximation  to  Ui  was  picked  so  that  it  had  no  negative 
value  in  the  domain.  The  approximation  to  U2  was  picked 
so  that  it  had  negative  and  positive  values  in  the  domain 
and  In  addition  it  was  orthogonal  to  Ui.  The  initial  ap¬ 
proximations  to  U3  was  picked  in  the  same  way  as  the  one 
for  U2  except  that  it  had  to  be  orthogonal  to  both  Ui  and 
U2.  It  may  be  difficult  to  make  a  suitable  choice  for  some 
of  the  higher  modes. 
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II .  RELAXATION  METHOD 


For  many  years  relaxation  techniques  have  been  used 
to  solve  differential  equations  with  and  without  the 
aid  of  computers.  They  are  basically  Iterative  proce¬ 
dures  in  which  a  new  approximation  is  obtained  from  a 
previous  approximation  and  its  residuals. 


A.  DERIVATION  OF  EQUATIONS 

In  this  section  a  typical  problem  is  posed  and  solved 
by  a  relaxation  method. 

Suppose  the  problem  to  be  solved  has  the  following 

form 


Z  +  Z 
xx  yy 


in  A, 


(2.1) 


where  the  unknown  function  Z  must  satisfy  the  differential 
equation  in  a  simply  connected  region  A  in  the  xy  plane, 
and  for  t>0.  In  addition  the  function  Z  is  required  to 
vanish  at  points  on  the  boundary  6A  of  the  region  A, 

Z(x,y,t)  =  0  on  6 A.  (2.2) 


A  typical  problem  is  that  of  a  vibrating  membrane.  The 
function  Z(x,y,t)  denotes  the  vertical  displacement  of 
the  membrane. 

Now  assume  that  the  displacement  has  a  representation 
of  the  form 


Z(x,y,t)  =  T(t)  U(x,y ) . 


(2.3) 
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When  Eq.  (2.3)  is  combined  with  Eq.  (2.1)  the  variables 
may  be  separated  ana  a  new  equation  is  obtained  of  the 
form 


(jit 

T“ 


+  U  ) 

yy. . 


u 


-x4. 


(2.1<) 


This  is  equivalent  to  two  equations 

T"  +  X2T  =  0  (2.5) 


and 


U 


xx 


*  \y  -  -x!u* 


(2.6) 


each  with  a  parameter  X.  Further  U  must  satisfy  the 
boundary  condition 

U(x,y)  *  0  on  6 A.  (2.7) 

This  is  a  typical  eigenvalue  and  eigenfunction  problem: 
to  find  the  values  X ,  or  X  ,  and  the  associated  functions 
U,  or  U^,  satisfying  Eq.  (2.6)  and  the  constraint  (2.7). 

The  values  Xn,  « i  a  X2  <  X3  <  are  called  the 

eigenvalues.  With  each  eigenvalue  is  an  eigenfunction 
U  .  In  this  problem  the  eigenfunctions  associated  with 
different  eigenvalues  are  orthogonal  on  the  region  A. 

Each  eigenfunction  is  also  called  a  mode. 

Consider  a  thin  elastic  membrane  of  a  particular  form, 
such  as  a  very  thin  uniform  sheet  of  rubber.  Assume  that 
the  membrane  is  made  fast  at  the  boundary,  while  it  is 
tightly  stretched  over  the  region  with  uniform  tension. 

Also  assume  that  damping  is  negligible.  Then  if  an  interior 
region  of  the  membrane  is  pushed  in  a  direction  perpendicular 
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to  the  plane  of  equilibrium,  it  becomes  distorted  into  a 
curved  surface.  The  resulting  area  can  be  computed  as 

S  -  /// 1  +  Ux*  +  Uyr  dxdy ,  (2.8) 

A 

Assume  that  U  and  U  are  very  small;  Eq.  (2.8)  becomes 
x  y 

approximately 

S  -  //( 1  +  *5UX*  +  *sUy2)  dxdy.  (2.9) 

A 

The  increase  in  the  area  of  the  membrane  due  to  the  distor¬ 
tion  is  therefore  approximated  by 

6S  »  //(l  +  %U  2  +  *sU  2 )  dxdy  -  //  dxdy  (2.10) 
x  y 

A  A 

-  h  /J(UX2  +  Uy2)  dxdy. 

A 

Hence  the  potential  energy  of  the  membrane  in  the  deflec¬ 
ted  position  is 

PE  =  v/2  //( U  2  +  U  2)  dxdy,  (2.11) 

a  y 

A 

where  v  is  the  tension,  assumed  to  be  constant  over  the 
region. 

Now  consider  any  particular  eigenfunction  or  mode. 

It  follows  from  the  solution  of  Eq.  (2.4)  that  the  deflec¬ 
tion  is  a  periodic  function  of  time  and  may  be  expressed 
in  the  form 

Z(x,y,t)  =  U(x,y)  sin  At,  (2.12) 

except  for  a  phase  shift.  Thus  Eq.  (2.11)  can  be  rewritten 
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as 


PE  *  h  JJ(U  2  +  U  2)  dxdy  sin2  Xt.  (2.13) 

A  Jr 

A 

The  maximum  value  of  the  potential  energy  is 

PEmax  *  v/2  //(Ux*  +  V5  dxdy*  (2.1*0 

A 

The  kinetic  energy  of  an  element  dm  ■  pdxdy  of  the  membrane 

Is 

*5pdxdy(Ut)2  -  *jpdxdy(U2x2  cos2xt),  (2.15) 

where  p  denotes  the  mass  per  unit  area  of  the  membrane. 
Therefore,  the  kinetic  energy  of  the  vibrating  system  is 

KE  ■  *sx2p  J/U2dxdy  cos2xt,  (2.16) 

A 

and  the  maximum  value  of  the  kinetic  energy  is 

KEjjjax  *  J/U2dxdy.  (2.17) 

A 

If  it  is  assumed  that  the  energy  is  constant  then,  the 
maximum  values  of  the  potential  and  kinetic  energy  are 
equal  for  individual  modes,  and  therefore 

*5X2p  //U2dxdy  =  v/2  //( Ux2  +  Uy2)  dxdy  (2.18) 

A  A 

or, 

v  //( U  2  +  U  2)  dxdy  (2.19) 

a  j 

X s  =  - - 

p  //U2dxdy 

A 
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The  first  eigenfunction  Ui  is  the  function  which  minimizes 
this  quotient  and  the  first  eigenvalue  Xj2  is  the  corres¬ 
ponding  value  for  the  quotient.  The  next  eigenfunction 
is  the  function  U2  which  minimizes  the  quotient  in  the 
space  of  functions  orthogonal  to  Ui,  and  X22  is  the  corres¬ 
ponding  value  of  the  quotient.  The  third  eigenfunction 
U$  minimizes  the  quotient  in  the  space  of  functions  ortho¬ 
gonal  to  Ui  and  U2,  etc.  The  set  U*,  U2,  ...  is  unique 
except  that  if  two  or  more  eigenfunctions  have  the  same 
eigenvalue,  they  may  be  replaced  by  linear  combinations 
of  themselves.  It  is  seen  that  Xi  <  X2  <  Xj  ... 

If  the  mode  U  is  known,  Eq.  (2.19)  can  be  used  to  obtain 
X2.  An  interesting  fact  is  that  a  rather  poor  approximation 
to  the  first  mode  Ui,  chosen  to  satisfy  the  appropriate 
boundary  condition,  will  yield  a  surprisingly  good  estimate 
for  X12.  This  result  is  apparently  due  to  Lord  Rayleigh, 
and  the  quotient  is  sometimes  called  Rayleigh* s  quotient. 


B.  COMPUTATIONAL  ROUTINE 

In  this  section  the  results  of  section  A  will  be  used 
to  develope  a  method  to  solve  Eq.  (2.6). 

If  Eq.  (2.6)  is  expressed  as  a  difference  equation 
then  it  may  be  written  as 


(UI+1,J  +Ui-l,J  -  2Ui,J> 


(Ui,j+l  +  Ui,j-1  '  2Ui,j} 


(2.20) 
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where  h  and  k  denote  the  x  and  y  mesh  size  respectively.  If 
the  mesh  sizes  sure  equal  then  Eq.  (2.20)  can  be  rewritten  as 


U 


(Ui+l,J  +  Ui,J+l  +  Ui-l,J  +  Ui,J-l)  C2-2D 


1,3 


4  - 


The  procedure  used  to  solve  the  partial  differential 
equation  was  to  pick  a  function  U°  which  satisfied  the  given 
boundary  conditions  as  a  first  approximation  to  U(x,y). 

This  function  need  not  be  a  very  good  approximation  to  U, 
and  in  fact  step  functions  were  sometimes  used. 

Prom  Eq.  (2.19)  an  approximation  to  X2  was  obtained 
by  assuming  that  the  approximation  U°  was  the  desired  func¬ 
tion  U.  With  this  first  approximation  to  X2  Eq.(2.21) 
was  used  to  obtain  an  improved  estimate  of  U.  The  form 
of  Eq.  (2.21)  used  to  obtain  the  vth  approximation  was 

V— 1  v-l  v  V 


(U 


U 


1,3 


*  Ul,3+l  *  Ul-l,3  *  Ui,J-l)  J  (2,22) 
”  (4  -  X2h2) 


in  which  the  v-lst  estimate  of  X  was  used.  By  alternating 
Eq.  (2.19)  and  Eq.  (2.22)  a  close  approximation  to  X2  and 
U  were  obtained.  The  solution  converged  to  the  smallest 
eigenvalue  Xi  and  the  associated  eigenvector  Ux . 


C.  HIGHER  ORDER  EIGENVALUES 

In  this  section  the  method  is  extended  to  find  the 
larger  eigenvalues  and  eigenfunctions. 

To  obtain  the  larger  eigenvalues  and  eigenvectors  the 
same  procedure  as  that  in  section  B  was  used  except  that 


t 


additional  equations  were  added  to  force  the  eigenfunctions 
to  be  orthogonal  to  those  already  found.  Define  for 
i  *  1,  2,...,  n  to  be  the  ith  eigenfunction,  associated 
with  X±.  All  higher  eigenfunctions  must  be  orthogonal  to 
the  lower  ones  obtained.  Orthogonalization  was  accomplished 
by  subtracting  out  multiples  of  the  lower  eigenfunctions 
already  obtained.  This  was  effected  by  expressing  U  as 


U 


i+1 


U 


new 


i+1 


old 


n 

-  I  c,u., 

j-1  J  J 


where 


J/U1+l  U1  <*xdy 

A  1  xold  J 

C  *  - 

J  //U*  dxdy 
A  J 


J-1, 2 


(2.23) 


(2.24) 


For  each  eigenfunction  desired  a  different  initial 
approximating  function  was  used.  The  method  gave  the  eigen¬ 
values  and  associated  eigenfunctions  in  numerically  increas¬ 
ing  order.  The  method  was  very  simple  and  effective  for 
lower  eigenvalues . 

The  fault  of  the  method  is  that  convergence  was  poor 
and  computation  times  were  large  if  the  number  of  intervals 
in  both  directions  was  large. 
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III.  DYNAMIC  PROGRAMMING  £8,10j 


The  method  of  dynamic  programming  has  recently  been 
applied  to  the  solution  of  partial  differential  equations, 
[10].  The  difference  equation  may  be  regarded  as  leading 
to  a  system  of  linear  equations  with  a  large  number  of 
unknowns.  The  method  effectively  reduces  the  number  of 
unknowns  involved  so  that  a  number  of  systems  are  to  be 
solved,  ‘each  one  of  much  lower  dimension.  In  one  case, 
for  example,  using  a  grid  with  N+l  intervals  in  the  region 
in  each  direction,  N  systems  each  with  N  unknowns  are  solved, 
rather  than  one  system  with  N  -  squared  unknowns. 

In  section  A,  the  method  is  applied  to  the  solution  of 
Poisson’s  equation  over  a  rectangle,  following  the  paper  of 
Angel  ClOl-  In  section  B,  the  method  is  extended  to  a  re¬ 
lated  eigenvalue  problem  by  combining  it  with  Stodola's 
method.  In  section  C,  the  method  Is  extended  to  the  solu¬ 
tion  of  the  biharmonic  equation,  by  applying  dynamic  pro¬ 
gramming  twice.  Finally,  in  section  D,  the  method  is 
applied  to  the  eigenvalue  program  involving  V4U  =  X2U,  again 
by  combining  the  method  with  Stodola's  method. 

A.  DYNAMIC  PROGRAMMING  FOR  PARTIAL  DIFFERENTIAL  EQUATIONS 
Consider  the  solution  of  Poisson's  equation 

Uxx  +  Uyy  “  f(x>y)  in  A,  (3.1) 

where  U  =  U(x,y)  is  subject  to  given  boundary  conditions 
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such  as 


(3.2) 


U(x,y)  -  g(x,y) 

It  may  be  noted  that  Eq.  (3.1)  is  the  Euler  equation 
associated  with  the  variational  problem 

j 

I(U)  -  Min  //( U  2  +  U  2  +  2fU)  dxdy  (3-3) 

U  A  *  y 

* 

where  the  function  U  is  chosen  from  the  class  of  functions 
with  first  partial  derivatives  belonging  to  L  over  A,  and 
satisfying  the  boundary  conditions  of  (3*2)  on  5 A. 

Let  the  region  A  be  discretized  by  choosing  n+1  and 
m+1  equally  spaced  points  in  the  independent  variables 
x  and  y  respectively.  Then  Eq.  (3.3)  may  be  rewritten  in 
a  discrete  version  with  equal  intervals  as 


KU) 


**  Min  l  l  f(U.  .  -  U±  ^-1)2 

U,  .  i=l  1,J 

l  >  J 

+  <ui,j  -  +  2fi,j  ui,jh2]> 


(3.4) 


where  (u  , (U.  ),  (U  , },  and  {U,  >  are  determined 

o,j  i,o  n,J  i,m 

from  the  boundary  conditions  of  Eq.  (3.2).  Now,  if  all  terms 
in  Eq.  (3*4)  involving  only  boundary  values  were  removed, 
while  not  affecting  the  solution,  a  more  convenient  form 
of  Eq.  (3.^)  is  obtained 


I(U)  = 


n  m 

i  i 

1=1  j=i 


Vj-i5 


(3-5) 


m-i  m-i 

+  3k  2h2fi,juij  +  3k  (ui,j  - 
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In  vector  notation  Eq.  (3.5)  may  be  rewritten  in  the  form 


I(U)  »  Min  I  (<QUr,Or>  +  <rRJUR>  +  S 


UR  i-1 


R*  R 


(3.6) 


+  R  -  ^R-l*  ^R  '  °R-1  ^  +  <  ^R*^R>^ 

In  this,  Ur  *  (UR  . . ,  UR  )T.  This  relation  now  defines 

1  m-1 

a  symmetric  matrix  Q,  vectors  rR,  and  Fr,  and  a  scalar 
SR  by 


2  i  *  J 

-  {q±  where  qA  ^  *  -  -1  |i-j|  «  1 

0  otherwise 


(3.7) 


rR  *  <rR  <>»  where  rR  4* 


fR  »  {2h2fRjJ} 

SR  *  UR,o  +  UR,m  * 


“2UR,o  J  '  1 

"2UR,m  J  *  -1 
0  otherwise 


The  notation  as  in  <rR,  UR>,  denotes  the  scalar  product  of 
the  two  corresponding  vectors.  The  matrix  Q  remains  constant 
while  rR  and  SR  are  functions  of  the  boundary  conditions 
only. 

Now  in  order  to  solve  Eq.  (3*6)  a  sequence  of  dynamic 
programming  problems  was  considered.  Let 


Fp  (7)  =  Min  l  (<QU,  ,  U,>  +  <r,  ,  U, >  (3.8) 

R  Un. . .U  ,i=R  1  1  1  1 


'R - n-1 


+  Si  +  <^'i  -  ^i-l»  "  ^i-l>  +  <^'i  * 
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where  tTR_1  *  7  and  tTn  is  given  by  Eq.  (3.2).  Now  Eq.  (3.3) 
can  be  rewritten  as 

FR(V)  ®  Min  (<5ffR.ffR>  +  <rR,  ffR>  +  SR  (3.9) 

^R 

+  <UR  -  ^,UR  -  V>  +  <fR  ,  U"R>  +  <$tfR+1>  tfR+1> 


+...  +  Sn  +  <Vn 


vs. 


n-11 


n 


-  Vl>  +  «n‘  V>- 


This  may  be  done  since  the  minimization  over  UR+1,..., 

tT  ,  can  be  commuted  with  the  minimization  over  tT„.  This 
n-x  w 

is  a  common  technique  of  dynamic  programming  based  upon  the 
principle  of  optimality,  £83.  This  allows  Eq.  (3-9)  to  be 
rewritten  as 


*R  L 


Fr(V)  »  Min [<QUr»  trR>  +  <rR>  trR>  +  SR  (3.10) 


+  <UR  -  7,  uR  -  V>  +  <fR>  uR>  +  FR+1(UR) 


'R 


R *  WR  R+l'  R‘ 


The  final  equation  is 


Fn(V)  *  <QUn>  u„>  +  <rn>  un>  .  Sn 


(3.11) 


+  <un  -  V-  Un  -  V>  -  <rn*  V • 


and  Un  is  known  from  the  boundary  conditions  Eq.  (3.2). 
Since  FR(V)  is  quadratic  in  V  Eq.  (3.11)  may  be  rewritten 
in  the  form 


Fr(V)  =  <ARV,  ?>  +  <Fr7>  +  CR.  (3-12) 

Now  by  substituting  from  Eq.  (3.12)  into  Eq.  (3.10)  and 
then  differentiating  with  respect  to  UR  an  expression  for 
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is  obtained  of  the  form 


°R  *  <r  +  5  +  Vi>rl 


’R+l  +  rR  +  fR 


.  ' (3.13) 


This  is  obtained  by  substituting  rel'atecjl  Eq,  (3-13,)  into 
Eq.  (3.10)  and  then  combining  Eq.  (3.10)  with  Eq.  (3.12$, 
the  various  quantities  in  Eq.  (3*13)  are  defined  by  these 
steps  as  ; 

^R 


=  T  -  (I  +  Q  + 

A  )“1 

rtR+l * 

(3.1*0 

j 

•  (T+  «  +  W*  <bR-l  *  rR  +  ,fR>. 

(3.15) 

! 

CR+1  +  SR  “ 

fd  .  a  .  A^r1 

'  (3.^6) 

<  Vl  4  rR,  bR+l  +  rR 


with  initial  values  determined  from  Eq.  (3.11) 'as' 


■  T- 


bN  *  -2trn* 


Cn  =  <(T  +  8)  (Tn>  V  +  <rn  -  V  +  Sn. 


(3.17) 

1 

(3.18) 

(3.19) 


The  matrix  (T  +  Q  +  ar+^)  is  nonsingular,  1  .  Thus  It  has 
inverses  which  may  be  computed  beforehand,  since  it  is: de¬ 
pendent  only  upon  the  type  of  operator. 

Due  to  the  fact  that  only  the  values, of  UR  are  desired 
the  quantities  CR  need  not  be  calculated. 

The  procedure  is  to!  calculate  the  quantities  in 
Eq.  (3*1*0  repeatedly  until  Aj  and  5*2  are  obtained.  '  Then 
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TX  is  found  from  Eq.’  (3.13)  with  V  ®  U  which  is  known 

•  f  o 

from  the  boundary  conditions .  Next,  Eq.  (3.13)  is  repeated¬ 
ly  solved  using  the  stored  values  of  ftR  and  B"R>  and  the 
last  value  of  tJR  as  V. 

'  I 

Thus,  the  problem  was  solved  by  n-1  inversions  of  sym- 

I  ; 

metric  matrices  of  order  m-2.  While  these  matrices  may  be 
large  thdre  are  efficient  computer  routines  available  to 
determine  the  inverses.  Once  the  inverses  are  found,  they 
may  be  stored  for  future  use,  since  they  are  based  only  on 
the  geometry  of  the  region  A.  Thus  for  several  problems 

f 

over  the  same  ^eg^.on  the  inverses  may  be  entered  into f the 

i  4 

program  as  data.  Also  they  have  the  property  that  if  less 
than  n+1  gild  'points  are  required  a  reduced  number  of  the 
matrices  may  be  used. 

B.  DYNAMIC  PROGRAMMING  FOR  EIGENVALUE  PROBLEMS 

In  the  first  section  of  this  chapter  it  was  seen  that 

: 

dynamic  programming  could  be  used  to  solve  partial  differ¬ 
ential  equations.  In  this  section,  the  program  is  modified 
to  solve  a  related  eigenvalue  problem  by  Stodola's  method. 

Consider'  the  problem  of  the  vibrating  membrane  consi¬ 
dered  in  Chapter  II.  The  differential  equation  for  the 
eigenvalues  and  eigenfunctions  is 

i  ■ 

:Uxx  +  Uyy "*  “X*U  '  ln  A,  (3.20) 

;  I 

where  U  s  U(x,y)  is  subject  to  the  boundary  condition 

U(x,y)  *  0  on  6A.  (3.21) 
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This  problem  was  solved  by  two  related  methods.  In 
both  Eq.  (3. 20) is  regarded  as  a  special  case  of  Eq.  (3*1) 
in  which 


f(x»y)  s  -X*U(x,y)  in  A  (3.22) 

* 

and 


f(x,y)  -  0 


on  <5A  (3.23) 


1 .  First  Routine 

The  difficulty  here  was  that  the  function  f  was 
known  only  on  the  boundary  and  \  was  unknown.  The  first 
step,  to  overcome  this  obstacle,  was  to  choose  a  function 

e 

U°Ix,y) which  satisfied  the  given  boundary  conditions.  Then 
an  estimate  of  x1  was  obtained,  say  (x°)2,  by  using  Ray- 
leigh's  formula.  Next  a  new  estimate  of  U,  say  UJ(x,y),  was 
obtained  using  dynamic  programming  to  solve  the  equation 

uix  +  uyy  “  -(x°)2  u°  s  (3.24) 


By  repeating  this  sequence  of  Rayleigh's  formula  and  dyna¬ 
mic  programming,  a  good  approximation  to  the  minimum  eigen¬ 
value  and  the  associated  eigenfunction  was  obtained.  The 
next  two  eigenvalues  and  vectors  were  also  obtained  by  the 
process  of  forcing  the  higher  eigenfunctions  to  be  orthogo¬ 
nal  to  ones  already  obtained  as  was  done  in  relaxation. 

The  inverse  matrices  used  in  the  routine  were  calculated 
in  determining  the  first  eigenvalue  and  eigenfunction; 
they  were  then  stored  so  that  it  was  not  necessary  to  recal¬ 
culate  them  for  subsequent  eigenvalues. 
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2.  Stodola's  Method 


The  second  form  is  more  like  the  usual  form  of  Sto¬ 
dola's  method  and  hence  a  brief  of  Stodola's  technique  is 
given  first. 

An  initial  function  V0  satisfying  the  boundary  con¬ 
ditions  is  selected.  It  may  be  considered  to  be  of  the 
following  form 

V0  «  aiUi  +  aaU2  +...  (3.25) 

where  ai  is  not  equal  to  zero.  For  convenience  VQ  was 
normalized;  the  L„  norm  of  V0,  | |V0| ( ,  is  defined  as 

MAX  |V0(x,y)|  -  l'|V0]jj  (3.26) 

A 

and  ||V0| |  is  set  equal  to  one.  Now  consider 


V0 


(3-27) 


and 


■*m . 


f-1 

m 

M 

aiU i  +  a  2 


2m 


Uz  +. . . 


(3.28) 


In  Eq.  (3.28)  it  can  be  seen  that  the  relative  size  of  the 
components  U2,  Uj,...  are  decreasing  by  a  factor  (Xj/X2)2, 
(Xi/X3)2*...  respectively.  After  a  few  applications  of  the 
operator,  L**1,  to  V0,  the  leading  term  will  dominate. 

Thus  the  functions  obtained  will  approximate  Vj,  except 
for  a  constant  factor,  and  the  ratio  of  successive  iterates 
will  approach  a  constant,  -1/x2. 


3 •  Second  Routine 


This  was  applied  as  follows.  Let  2^  be  defined  by 
the  equation  LZm  -  V  50  that 


Zm  -  L-1V  . . 
m  m-i 


(3.29) 


This  equation  was  solved  by  dynamic  programming.  It  was 
found  convenient  to  normalize  each  Iteration.  Let 


(3.30) 


Then  the  approximation  can  be  made 


X;2  * 


m' 


(3.31) 


The  sequence  of  functions  V0,  V i,...  converges  to  Ui. 

The  process  was  terminated  when  successive  approx¬ 
imations  were  sufficiently  close  together,  say, 

||Vm(x,y)  -  vm_1(x,y) | |<e,  (3*32) 

where  e  is  a  preassigned  small  positive  number.  The  result¬ 
ing  function  Vm  is  an  approximation  to  Uj,  and  Xj  is  ob¬ 
tained  from  Eq.  (3.31)* 


C.  DYNAMIC  PROGRAMMING  FOR  A  HIGHER-ORDER  OPERATOR 

In  section  A  and  B  of  this  chapter  it  was  shown  how 
dynamic  programming  could  be  used  to  solve  Poisson’s  equa¬ 
tion  and  the  vibrating  membrane  program,  respectively. 

In  this  section,  by  another  modification  to  the  routine, 
the  biharmonic  equation  can  be  solved. 


2H 


Assume  the  problem  to  be  solved  is  of  the  form 


?*(U)  *  f (x,y)  in  A,  (3.33) 

where  U  «  U(x,y)  is  subject  to  the  constraints 

U  «  g(x,y)  on  <$A  (3.3*0 

32U/3n2  -  P(x,y )  on  6A. 

Clearly  Eq.  (3.33)  may  be  rewritten  in  the  form 

V1  4(x,  y)  «  f(x,y),  (3.35) 

where 


V2  U(x,y)  =  *(x,y).  (3.36) 

Now,  since  U  and  32U/3n2  are  both  known  as  the  boundary, 
4(x,y)  may  be  approximated  on  the  boundary.  On  a  rectangle, 
for  example,  the  following  relations  determine  $  on  the 
boundary 

4(0, y)  *  -P(0,y)  +  Uyy  (0,y )  (3.37) 

4(n,y)  =  p(n,y)  +  nyy  (n,y) 

4(x,0)  ^  -I ( x ,0 )  +  Uxx  ( x , 0 ) 

4 (x,m)  =  P(x,m)  +  Uxx  (x,m). 


The  usual  finite  difference  scheme  may  be  used  to  approxi¬ 
mate  U  and  U  and  thus  the  relations  of  Eq.  (3.37)  may 
ax  yy 

be  approximated  by 
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♦  .  *»  p  ,  + 

nj  nj 


(3.38) 


for  i  "  X  , . *  *  ,  n-1  j  J  X  , « . « ^  m-3. 


and  <$n  ,  4>n  « »  „  are  known  from  the  boundary  con- 

ditions.  Thus  <t(x,y)  is  now  approximated  on  the  boundary. 
First  Eq.  (3.35)  is  solved  by  dynamic  programming  to  obtain 
the  function  $  in  the  region.  Then  Eq.  (3.36)  is  solved 
by  dynamic  programming  to  obtain  the  desired  solution. 


D.  DYNAMIC  PROGRAMMING  FOR  HIGHER  ORDER  EIGENVALUE  PROBLEMS 
This  method  may  be  extended  to  the  corresponding  eigen¬ 
value  and  eigenfunction  problem  much  as  before.  One  such 
physical  problem  is  that  of  a  vibrating  uniform  plate  with 
hinged  edges . 

Assume  that  the  differential  equation  has  the  form 

V4U  *  X2U,  in  A,  (3.39) 

and  the  boundary  conditions  are 

U  =  —  =0  on  6A  (3*40) 

3n2 

The  technique  of  sections  B  and  C  may  be  applied  in  two 
steps  to  the  solution  of  the  problem.  Let  4>  be  defined  as 
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«  =  v2u. 


(3.41) 


Then  the  two  equations  to  be  solved  are 

a  x2Un  (3.42) 

and 

y2yn+l  s  (3.A3) 

subject  to  the  boundary  conditions,  Eq.  (3.^0) • 

In  order  to  start  the  routine  an  initial  estimate  for 
the  function  U(x,y),  say  (U±J>,  is  made.  The  value  of  X 
is  estimated  as  was  done  in  section  B.  Then  Eq.  (3^2) 
and  Eq.  (3.  *»3)  are  solved  by  dynamic  programming  to  get  the 
next  approximation  {U^J.  The  same  criterion  for  stopping 
is  used  as  that  in  section  B. 


IV.  RAYLEZOH-KITZ  METHOD  £^*5,6,9] 


The  Rayleigh-Ritz  method  has  been  used  for  many  years 
to  obtain  approximations  to  the  solution  of  partial  differ¬ 
ential  equations  and  eigenfunction  problems.  In  this  method 
the  problem  is  posed  as  a  minimization  problem,  say  invol¬ 
ving  an  integral.  Then  some  linearly  Independent  functions 
which  satisfy  the  boundary  conditions  are  chosen.  The  solu¬ 
tion  is  approximated  by  a  linear  combination  of  these. 
Finally  the  coefficients  in  the  approximation  are  chosen 
so  as  to  effect  minimization.  This  leads  to  an  eigenvalue 
problem  involving  symmetric  matrices. 

The  functions  chosen  may  be,  for  example,  polynomials 
of  low  degree,  or  trigonometric  functions.  Assume  that 
homogeneous  boundary  conditions  are  given.  Let  $k  =  $k(x,y) 
be  n  functions  which  satisfy  these  and  approximate  the  solu¬ 
tion  U  by 

U(x,y)  =  l  Ck*k.  (ii.l) 

k  1 

A.  PARTIAL  DIFFERENTIAL  EQUATIONS 
In  this  section  the  solution  of 

Uxx  +  Uyy  =  f(x’y)  in  A’  (1*3) 

is  again  considered.  It  is  the  Euler  equation  associated 
with  minimizing  the  integral 
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(4.2) 


J/(U  2  +  U  2  +  2fU)  dxdy. 
A  x  y 


Substituting  from  Eq.  (4.1)  into  Eq.  (4.2)  and  integration 
results  in  a  function  which  may  be  written 

I  -  I(C1,...,  Cn),  (4.3) 


for  functions  of  the  form  (4.1).  The  minimization  of  this 
function  and  an  approximate  solution  to  Eq.  (1.3)  is  thus 
obtained  by  solving 
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0, 


k  a  1,  2, . . . ,  n. 


(4.4) 


The  effectiveness  of  the  procedure  of  course  depends  upon 
the  choice  of  the  approximating  functions 


B.  EIGENVALUE  PROBLEMS 

The  method  is  also  applicable  to  eigenvalue  problems. 
Consider  again  the  problem  in  Chapter  II  of  finding  eigen¬ 
values  and  eigenfunctions  for  the  equation 

Ux*  +  Uyy  “  ln  A  <2-6> 

subject  to  the  boundary  condition 

U(x,y)  =0  on  &A.  (2.7) 

2 

In  many  problems  Xi,  the  lowest  eigenvalue,  is  the  minimum 
of  the  ratio  of  two  integrals.  This  fact  was  shown  in 
Chapter  II  and  Eq.  (2.19). 

If  an  approximation  to  U  is  chosen  as  in  section  A,  since 
each  function  satisfies  the  linear  homogeneous  boundary 
condition,  the  sum  shown  in  Eq.  (4.1)  satisfies  it.  If 
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this  sum  is  substituted  into  Eq.  (2.19)  the  resulting  values 

define  an  upper  bound  for  X,  for  all  choices  of  the  C*s. 

Further  also,  the  C's  may  be  chosen  to  give  a  least  upper 

bound  over  the  subspace  spanned  by  $x,  4a,.*.,  $n. 

If  that  approximation  for  U  is  substituted  into  Eq.(2.19) 

the  numerator  and  denominator  become  quadratic  forms  in 

(C  C  )  =  C.  The  numerator  has  the  form 

1  n 

I  ?  »1,  qc,  *  ?  <*.5> 

i-1  J«1  1  J 

where 

a.i  «  //(^  +  4.  ♦  ,  )  dxdy,  (4.6) 

J  A  x  Jx  xy  Jy 

and  the  denominator  has  the  form 


l  l  b  C  C  .  V 
i«l  J«1  1  J 


(4.7) 


where 


ij 


//* i 

A  1  J 


dxdy . 


(4.8) 


These  define  the  matrices  A  and  B.  The  first  eigenvalue 
X x  minimizes  the  quotient  in  Eq.  (2.19)  and  hence  the  mini¬ 
mum  value  of 

^fcA  C/C^B  C 

defines  the  minimum  value  of  the  quotient  in  the  subspace 
spanned  by  $  x , .  . . ,  4>n* 
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To  find  this  is  equivalent  to  minimizing  the  quadratic 
form 

Xi*  -  Min  C*!  C  (4.10) 

S’ 

subject  to  the  constraint  that  the  second  quadratic  form 
assumes  the  value  one 

3*5  £-1.  (4.  11) 

The  problem  may  be  solved  in  two  steps.  First  find  the  eigen 

values  and  eigenvectors  of  B.  Let  the  eigenvalues  of  B 
be  i*l,  2,...,  n,  and  the  associated  normalized  eigen¬ 

vectors  V^.  Let  T  be  the  matrix 

f  -  (Vi,  V*,...,  Vn)  diag(l/«i ,  l/wn).  (4.12) 
Then  the  transformation 

C  -  T  D  (4.13) 

reduces  the  constraint  (4.11)  to  the  form 

^5  S’  *  -1.  (4.14) 

condition  (4.10)  becomes 

Xi2  =  Min  D,  (4.15) 

n 


where 

B  -  T,  (D. 16) 

subject  to  the  constraint  (4.14). 

Hence  the  value  of  Xi2  is  the  smallest  eigenvalue  of 
the  matrix  E. 
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Let  Ci  be  the  associated  normalized  eigenvector  of  E-. 
Then  the  corresponding  minimizing  function  has  coefficients 
determined  from  Eq.  (4*13) 

Ci  *  ¥  Di.  (4.17) 

This  value  of  Xi2is  an  upper  bound  for  the  first  eigen- 
value,  the  least  upper  bound  in  the  space  of  functions 
spanned  by  $1,...,  <5^.  In  a  similar  way  the  second  eigen¬ 
value  of  E  furnishes  an  upper  bound  to  the  second  eigenfunc 
tion,  etc. 

For  simple  problems,  at  least,  it  seems  easy  to  choose 
the  4*  *  s  so  that  a  good  bound  for  the  lowest  eigenvalue  re¬ 
sults.  It  is  not  clear  how  to  make  good  choices  to  get 
good  estimates  for  the  higher  eigenvalues. 


* 
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V.  QALERKIN’S  METHOD  [H,5,6j 


Sometimes  the  solution  to  boundary-value  problems  in¬ 
volving  partial  differential  equations  can  be  obtained  by 
forming  an  associated  functional  in  the  form  of  a  definite 
integral  which  is  to  be  made  stationary.  For  a  problem  of 
this  type,  the  Rayleigh-Ritz  method  is  often  an  effective 
procedure  for  the  determination  of  an  approximate  solution. 
However,  in  many  instances  it  is  difficult  to  find  this 
functional.  In  such  cases,  the  Galerkin  method  is  often 
effective. 

A.  PARTIAL  DIFFERENTIAL  EQUATIONS 

Consider  the  linear  homogeneous  boundary  value  problem 
of  the  form 

LU(x,y)  =  f (x,y) ,  (5.1) 

subject  to  linear  homogeneous  boundary  conditions.  The 
symbol  L  stands  for  a  linear  differential  operator,  such 
as  v2. 

Suppose  that  an  approximate  solution  is  taken  in  the 

form 

U  (x,y)  =  l  C.  4>.  (x ,y ) ,  (5.2) 

n  k-1 

where  the  coefficients  C^,...,  Cn  are  constants,  to  be  de¬ 
termined,  and,  as  in  the  previous  chapter,  the  functions 
4>i»...,  $n  are  picked  to  satisfy  the  homogeneous  boundary 
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conditions.  The  coefficients  are  dependent  upon  'the  number 
of  functions  picked  and  therefore*  must  be  recomputed  if 
a  larger  number  of  functions  are  later  chosen.  While  the 
function  Un  satisfies  the  boundary  conditions*  it  will1  not, 
in  general,  satisfy  Eq.  (5.1)i  Thus  there1  is  a  residual,  1 

LUn(x,y)  -  f(x,y)  “  Rn(x,y),  (5.3)  . 


which  can  be  viewed  as  an  error  or  penalty  function.  Assume 

I 

the  solution  for  U(x,y)  can  be  expressed  by  an  infinite 
complete  series  of  these  linearly  independent  functions 
in  the  form  1 


U(x,y)  *  I  ck*k(x»y) 
k=l  K  K 


<5.M) 


Now  Un(x,y)  in  Eq.  (5.2)  represents  a  sequenc'e  of  partial 
sums  which  approximate  U(x,y).  Now  if  the  condition  is 

i 

imposed  that  L(U  )  -  f  be  orthogonal  to  each  function 
*j_(x,y)  on  the  demaln  A,  the  following  set  of  equations 
are  obtained  (55) 

//( L(U  )  -  f)  t.  dxdy  =  0 ,  for  k  »  1',  2 , . . . ,  n 

a  n  k.  . 


or  by  use  of  Eq.  (5.3) 


//Rn(x*y)  $k(x,y)  dxdy  s  0,  for  k  =  1,  2 
A 


•  (5.6) 

. n. 


If  Eq.  (5.6)  is  to  hold  as  n  -*■  ®  it  follows  that 


lim  R  -  0 
M  n 
n-*-~ 


(5.7) 


3^ 


I 


I 


I 


Let  R(x,y)  be  an  arbitrary  function  satisfying  the 

i 

homogeneous  boundary  conditions.  Since  for  k  =  1,  2,... 

'  t  *  * 

forms  a  complete  set  of  functions,  constants  a,,...,  a  . . . . 

•  *  -L  n 

i 

can  be  found  such  that 

f 

*  J  GO  • 

R(x,y)  a  l  ak4>k(x,y).  (5.8) 

i  k-1,  K 

i 

Thus, 

I  ‘ 

//  lim  R  ( x ,  y )  n(x,y)  dxdy  *  0,  (5.9) 

A  n+“.  n 

i  1  •. 

for  any  arbitrary  n(x,y)!,  so  that  by  the  fundamental  lemma 
of  variational  calculus  Eq.  (5.7)  is  true  "almost  every¬ 
where".  Suppose  further  that  L(Un)  L(U)5then  Eq.  (5.1) 
holds,  by  the  use  of  Eq.  (5.7). 

Qalerkin's  method  requires  that  the  error  function 
Rn(x,y)  be  orthogonal  to  each  of  the  functions  <t>k,  or  that 
Eq.  (5.5)  and  Eq;  (5*6)  hold.  Now  by  substituting  Eq.  (5*2) 
into  Eq.  (5.5)  an  integral  is  obtained  of  the  form 

I 

f 

,  //  L(  l  0k$k)  -f  ^  dxdy  =  0,  (5.10) 

A  k=l 

i  i 

i  «  1 , '  2 , . . . ,  n . 

I 

This  results  in  a  system  of  n  linear  algebraic  equations 
in  n  unknowns  C^,...  C  ,  Furthermore,  the  system  is  in¬ 
homogeneous  unless  the  function  f(x,y)  is  orthogonal  to 
each  4>k(x,y);  Eq.  (5.10)  may  be  rewritten 


t 
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(5.11) 


?  C.  // L(*.)  dxdy  *  JJf*k  dxdy, 
J-l  J  A  J  K  A  * 


it  1  ,  2  ,  .  .  .  ,  ft. 

In  matrix  notation  this  becomes 


B  C  a  F, 


where 


C  - 

B  - 

F  = 


(Cl  C2  ...Cn)T 

£  to  ^ }  X  =  X)  2 , . . . ,  ti  ^ 
J  —  1>  2  j » . • |  n 


(fi***.j  fn) 


T 


and 


bij  “  //L,(4>i^  *j  dxdy 
A 

f.  a  dxdy 

1  A  1 


(5.12) 

(5.13) 

(5.14) 

(5.15) 

(5.16) 

(5.17) 


The  constants  are  obtained  by  solving  the  system  Hq. 

Eq.  (5.12). 

Collocation,  a  convenient  variation.  One  way  to  get 
an  approximate  solution  for  the  C’s  is  the  following. 

Choose  n  points  of  the  region  k.  Evaluate  the  terms  in 
the  integral  of  Eq.  (5-10)  at  these  points.  This  yields  n 
equations  for  the  unknowns  C  .  This  method,  called  colloca¬ 
tion,  is  not  generally  so  accurate  but  it  is  quicker  than 
carrying  out  the  integrations  to  define  the  coefficients 
Eq.  (5.16)  and  Eq.  (5.17) ■ 
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Qalerkins  method  Is  also  useful  to  obtain  a  direct 
solution  of  variational  problems.  Assume  It  Is  desired  to 
minimize  a  functional  of  the  form 


I(V)  -  //F(x,y  V,V  ,V  )  dxdy . 
A  *  y 


(5.18) 


It  Is  desired  to  find  an  extreme  value  for  the  functional 
subject  to  the  condition  that  V(x,y)  is  prescribed  on  the 
boundary  <5A  of  the  domain  A.  It  is  known  that  if  the  exist¬ 
ence  of  an  extremizing  function  U(x,y)  is  assumed  and  that 
the  f'-iction  F  possesses  continuous  derivatives  of  the 
second  order  with  respect  to  its  arguements,  there  results 
the  condition 


//b(x,y) 

A 


PU  ‘ 


3  p 
3“  FU 


x 


3“  FU 

y  y, 


dxdy  =  0 


(5.19) 


for  an  arbitrary  function  n(x,y)  which  has  piecewise  con¬ 
tinuous  derivatives  in  A  and  which  vanishes  on  <$A.  This 
is  derived  by  considering  a  function  of  the  form 


V(x,y)  =  U(x,y)  +  en(x,y), 


(5.20) 


and  differentiating  with  respect  to  e.  If  n(x,y)  is  other¬ 
wise  arbitrary,  then  one  form  of  the  fundamental  lemma  of 
variational  calculus  requires  that 


Suppose  now  that  n(x,y)  is  the  kth  function  <t^(x,y)  and 
that  an  orthogonality  requirement  is  imposed  as 
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y  VPu  "  I?  Pux  -  w  PV  *  °* 


(5.22) 


for  k  =  1,  2,. ...  n. 

Finally,  assume  that  the  function  U(x,y)  which  effects  the 
minimization  can  be  represented  satisfactorily  by  a  finite 
series 


n 

un(x,y)  *  l  ck4R(x,y).  (5.23) 

lc— ■  1 

Eq.  (5.22)  then  defines  a  system  of  n  algebraic  equation 

to  be  solved  for  the  n  unknowns  c,  ....  c  .  Thus  Galerkin's 

1  n 

method  is  applicable  in  the  solution  of  variational  problems. 
However,  much  of  its  value  lies  in  the  fact  that  it  is  not 
necessarily  connected  with  a  variational  procedure. 


B.  EIGENVALUE  PROBLEMS 

Suppose  that  it  is  desired  to  solve  an  eigenvalue  prob¬ 
lem  of  the  form 

LU  »  XU(x,y)  in  A,  (5-24) 

and 


U  =  0 


on  6A .  (5.25) 


Assume  as  in  the  first  section  that  the  function  U(x,y) 
can  be  approximated  in  the  form  of  Eq.  (5.23).  The  problem 
is  solved  by  considering  equations  of  the  form 


J/(L<V  -  XUn 


dxdy  s  0 


(5.26) 


for  J  -  1,  2,...,  n. 
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Now  by  substituting  for  Un 
tained  of  the  form 

n 


I 

J-l 


(ai,r  XYi,j)Cj 

for  i  =  1,  2,.. 


the  functions 


*  0 

•  »  n» 


a  system  is  ob- 


(5.27) 


where 


ai  J  =  J/lC^)  ♦j  dxdy  (5.28) 

A 

Yi,j  “  /J*i*j  dxdF  (5.29) 

^  A 


This  has  nontrivial  solutions  for  the  Cfs  if,  and  only  if, 

A  »  XYiJjls=  °»  (5.30) 

where  A  is  the  determinant  of  the  matrix.  Thus  the  values 
of  X^  may  be  found  by  solving  the  characteristic  equation 
(5.30).  For  each  eigenvalue  Xi  there  corresponds  a  system 
of  equations  (5.28)  to  be  solved  for  the  eigenvector 


J1  (ak,j  “  XiYk>j)  Ck  =  0 


(5. 3D 


for  J  s  1,  2,...,  n. 

Now  for  each  eigenvalue  X.j.  this  system  of  n  homogeneous 
equations  may  be  solved  to  give  the  values  where 

this  coefficient  corresponding  to  the  ith  eigenvalue. 

Thus  the  eigenvalues  are  obtained,  with  their  correspond¬ 
ing  eigenvectors.  The  functions  should  be  chosen  to 
reflect  whatever  characteristics  the  solution  is  felt  to 
have. 
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VI.  DISCUSSION  OF  VARIOUS  METHODS 


FOR  SOLVING  EIGENVALUE  PROBLEMS 


In  Chapters  II,  III  and  IV,  three  methods  were  developed 
for  finding  the  first  three  eigenvalues  and  eigenfunctions. 
All  three  gave  satisfactory  values  for  the  eigenvalues  and 
eigenfunctions,  and  satisfactory  times  for  the  rather  simple 
problems  considered  here.  These  methods  were  used  to  obtain 
solutions  of  the  following  two  problems: 

Ux*  -  uyy  -  -  x‘v  in  A- 

subject  to  the  constraint 

U(x,y)  “  0  on  6A;  (6.2) 

and  ■ 


V^U  *  X2U  in  A,  (6.3) 

subject  to  the  constraints 

U(x,y)  s  — -  =0  on  6A.  (6.11) 

3n2 


In  this  chapter  the  computation  and  numerical  results  are 
discussed  and  compared. 

In  all  of  the  methods  a  set  of  functions  was  needed. 

In  the  Rayleigh-Rltz  method  these  were  the  basis  for  tht? 
approximating  functions;  In  the  other  two  methods  they 
were  the  initial  estimates  of  the  functions.  The  ones  usu¬ 
ally  chosen  were 
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Ui°  -  (X  -  x2)(y  -  y2), 


(6.5) 


Ua°  a  (x  -  x2)(y  -  y2)(»s  -  x ) ,  (6.6) 

U*°  *  (x  -  x2 ) (y  -  y2)(»s  -  y).  (6.7) 

Step  functions  were  also  used  as  first  estimates  in  the 

iterative  methods;  these  increased  the  number  of  iterations 
required  some,  but  not  much,  particularly  for  the  higher 
modes . 

Usually  the  range  of  the  independent  variable  was  di¬ 
vided  into  ten  equal  sub-intervals ,  in  going  to  a  differ¬ 
ence  equation.  In  the  Rayleigh-Ritz  method  this  did  not 
yield  sufficient  accuracy  and  it  was  found  necessary  to  go 
to  forty  Intervals.  Twenty-five  intervals  were  also  used 
in  some  computations;  the  intermediate  number  was  chosen 
because  of  storage  requirements  in  dynamic  programming. 

In  the  iterative  procedures  some  convergence  or  stop 
criterion  was  needed.  When  a  relaxation  procedure  was  used 
together  with  Rayleigh's  formula  for  estimating  the  eigen¬ 
value,  computation  was  terminated  whenever  the  eigenvalue 
did  not  decrease  by  at  least  0.002.  In  Stodola's  method, 
the  routine  was  terminated  whenever  the  norm  of  the  change 
in  the  function  U(x,y)  was  less  than  0.002. 

The  relaxation  method  required  the  least  time  for  this 
simple  problem.  Most  of  the  time  required  for  dynamic  pro¬ 
gramming  was  spent  in  Inverting  the  matrices.  Dynamic  pro¬ 
gramming  yielded  a  very  accurate  approximation  to  the 
eigenfunction. 


Hi 


For  some  reason  it  was  necessary  to  use  forty  intervals 
to  get  the  desired  accuracy  with  the  Rayleigh-Ritz  method. 
When  ten  intervals  were  used,  the  eigenvalues  of  the  matrix 
I?  were  significantly  too  small,  apparently  due  to  errors 
in  the  integration  and  transformation  routines.  There  were 
at  least  two  other  disadvantages  of  the  Rayleigh-Ritz  method. 
First  it  was  tedious  to  program.  Second,  there  may  be  some 
difficulties  in  choosing  the  functions  particularly 
if  some  of  the  higher  modes  are  desired.  The  results  of  a 
set  of  computations  are  given  in  Computer  Output  5  ,  in¬ 
cluding  the  three  eigenvalues,  the  associated  coefficients 
and  the  values  of  the  corresponding  functions  at  various 
points.  The  routine  is  shown  in  Computer  Program  5.  The 
necessary  matrix  transformations  and  solutions  for  the 
eigenvalues  were  carried  out  using  programs  TRED2  and  TGL2 
respectively,  [7]. 

The  simple  relaxation  method  of  Chapter  II  had  the  ad¬ 
vantage  of  being  the  simplest  to  program  and  to  run.  It 
had  the  disadvantages  that  terminal  convergence  was  slower 
than  in  the  method  of  dynamic  programming  and  if  a  large 
number  of  points  were  involved  the  computing  time  increased 
greatly.  The  results  ** **  a  set  of  computations  are  given 
in  Computer  Output  1.  The  routine  is  shown  in  Computer 
Program  1. 

Dynamic  programming  had  the  advantage  of  yielding  very 
accurate  values  in  a  small  number  of  iterations.  The  dis¬ 
advantages  were  that  it  was  relatively  difficult  to  program 
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and  th~t  it  took  quite  a  bit  of  time  to  generate  the  inver¬ 
ses  of  the  matrices  required.  Computer  Output  7  shows  the 
results  by  this  method  and  the  program  used  is  in  Computer 
Program  6. 

The  values  of  the  eigenvalues  obtained  by  the  three 
methods  are  compared  in  Table  X.  In  Table  I  are  the  eigen¬ 
values  for  the  first  differential  equation,  the  difference 
equation  with  ten  intervals,  together  with  the  results  of 
he  computations,  the  number  in  parenthesis  by  an  entry 
indicates  the  number  of  intervals  used  in  the  computation. 
The  eigenvalues  for  forty  intervals  are  intermediate  between 
those  for  ten  and  those  for  the  differential  equation. 


Differential 

Equation 

Difference 
Equation (10) 

Rayleigh 

Relaxation 

Dynamic 

Programming 

4.44289 

4.  42463 

4.4110  (25) 

4.42486(10) 

4.42442(10) 

4.44751(40) 

4.44611(40) 

4.43753(25) 

7.02482 

6.92714 

7.23029(40) 

6.92736(10) 

6.92717(10) 

7 . 02482 

6.92714 

7.24707(40) 

6.92736(10) 

6.92717(10) 

Table  I. 

Comparison  of  First  Three  Eigenvalues  for  Eq.  (1.1) 


The  three  methods  gave  close  agreement  for  the  first 
eigenvalue  but  tended  to  differ  on  the  next  two. 

For  the  differential  equation  of  higher  order,  only 
relaxation  and  dynamic  programming  were  compared;  these 
comparisons  are  shown  in  Table  II  using  ten  intervals. 


Differential 

Equation 

Difference 

Equation(lO) 

Relaxation 

Dynamic 

Programming 

19.739227 

19.577361 

19.60767 

19.57777 

49. 348040 

47.985220 

48.00555 

47.98473 

49.348040 

47.985220 

48.02386 

47.98474 

Table  II. 

Comparison  of  First  Three  Eigenvalues  for  Eq.  (1.3) 


Computing  times  were  similar,  around  twelve  seconds 
for  each.  Computer  Output  2  shows  the  results  for  the  re¬ 
location  method,  and  the  program  is  Computer  Program  2. 
Computer  Output  9  shows  the  results  for  dynamic  programming 
and  the  program  is  Computer  Program  6 . 


VII.  DISCUSSION  OP  THE  SOLUTION  OP  A 
PARTIAL  DIFFERENTIAL  EQUATION  BY  VARIOUS  METHODS 

In  Chapters  III,  IV,  and  V,  three  methods  were  developed 
for  solving  partial  differential  equations.  These  methods 
were  used  to  obtain  solutions  of  the  following  problems : 

Uxx  +  Uyy  *  +  y*  -  x  -  y)  in  A,  (7.1) 

subject  to  the  constraint 

U(x,y)  -  0  on  6A;  (7.2) 

and 

V*U  *  8  in  A,  (7.3) 

subject  to  the  constraints 

U(x,y)  -0  on  6A 

=  2(y  -y2)  for  x  ■  0  (7.4) 

3n2 

=  -2  (y  -y2 )  for  x  =  1 

3n2 

— -  =  2 ( x  -x2)  for  y  =  0 

3n2 

— -  =  -2 (x  -x2)  for  y  =  0 

3n2 

In  this  chapter  the  computation  and  numerical  results  are 
discussed. 


In  Galerkin's  method*  an  approximation  for  the  function 
U(x,y)  satisfying  the  constraint  Eq.  (7.2)  was  made  by  choos¬ 
ing  suitable  functions  4>i(x*y),  $2(x,y)  and  *»(x*y)  to  satis¬ 
fy  the  constraint  in  Eq.  (7.2).  The  approximation  for  U(x,y) 
was 

U(x*y)  =  Ci$i(x,y)  +  C2$2(x,y)  +  C3<t>9(x*y)  (7.5) 

*  (x  -x2 ) (y  -y2)(Cj  +yxC2  +x2y2C$). 

The  values  of  Ci,  C2  and  C3  were  obtained  by  techniques 
described  in  Chapter  V.  Computer  Output  3  shows  the  values 
of  Ci,  C2,  Cj  and  the  function  U(x*y)  at  various  points. 

In  this  method  the  region  A  was  sub-divided  into  ten  equal 
sub-intervals  for  each  independent  variable  for  the  inte¬ 
gration  routine.  A  second  approximation  for  U(x,y)  was  made 
for  this  method  as 

U(x  *y )  =  Cx4i  +  02^2  +  C3<l>3  (7*6) 

*  (x  -x2)(y  -y2)(xCj  +  yC2  +  xyC3). 

The  purpose  of  this  was  as  follows.  There  generally  is  some 
skill  and  art  involved  in  choosing  the  functions  G>n 

well.  In  fact  in  the  first  choice  is  actually  the  desired 
function.  The  second  set  of  functions  was  chosen  so  as  to 
get  some  feel  for  the  consequences  of  a  poor  choice  of  the 
$^'s.  The  values  of  C*,  C2,  C3  and  the  function  at  various 
points  is  shown  in  Computer  Output  10. 

The  numerical  solution  of  Eq.  (7.1)  by  dynamic  program¬ 
ming  with  the  constraint  of  (7.2)  yielded  the  values  in 


Computer  Output  6.  Again  the  region  was  divided  as  was 
done  In  Galerkin’s  method. 

In  the  Rayleigh-Ritz  method,  an  approximation  for  the 
solution,  U(x,y),  satisfying  Eq.  (7.1)  was  made  by  choosing 
suitable  functions  $i(x,y),  #2(x,y)  and  *s(x,y)  to  satisfy 
the  constraint  (7.2).  The  approximation  for  U(x,y)  was 

U(x,y)  -  C i $ i  +  C2*2  +  Cj*3  (7.7) 

-  (x  -  x2)(y  -  y2)(Ci  +  xyC2  +  x2y2C3). 

The  values  of  Ci,  C2,  and  Cs  were  obtained  by  techniques 
described  in  Chapter  IV.  Computer  Output  4  shows  the  val¬ 
ues  of  Ci,  C2,  Ci  and  the  values  obtained  for  U(x,y)  at 
various  points.  In  this  method  it  was  found  necessary  to 
sub-divide  the  region  into  forty  equal  sub-intervals  for 
each  independent  variable  in  order  tc  get  satisfactory 
accuracy. 

The  best  approximation  to  U(x,y)  was  obtained  by  Galer- 
kin's  method  using  Eq.  (7*5),  where  #i  was  the  desired 
function.  There  was  no  error  and  the  method  was  able  to 
detect  that  this  was  the  case.  However,  when  Eq.  (7.6) 
was  used  the  maximum  error  was  eight  thousandths.  The  time 
required  to  solve  the  problem  by  this  method  was  0.57 
seconds . 

The  problem  was  solved  by  dynamic  programming  in  three 
seconds  with  accuracy  to  six  digits.  However,  over  half 
of  the  computer  time  was  spent  obtaining  inverses,  which 
could  have  been  fed  as  data  from  solving  Eq.  (6.1)  in  this 
particular  case. 
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The  Rayleigh-Ritz  method  required  Jyst  under  ten  seconds 
and  had  a  maximum  error  of  two  thousandths. 

i 

Because  of  the  time  required  and  the  accuracy  bbtained 
by  Rayleigh-Ritz ,  only  Galerkin’s  method  and  dynamic  pro¬ 
gramming  were  used  to'  solve  Eq  L  (7*3)..  : 

Galerkln's  method  obtained  the  same  values  and  accuracy  ' 
in  the  solution  of  Eq.  (7.3)  as  it  did  in  the  solution  to 
Eq.  (7.1)  and  took  the  same  time.  ’ 

Dynamic  programming  obtained  the  same  degree  of  accuracy 

I 

in  solving  Eq.  (7*3)  as  it  did  in  solving  Eq.  (7.1)  and 
took  three  seconds.  The  results  are  shown  in  Computer 
Output  8  and  the  program  is  Computer  Program  7. 

i  ■  ! 

l 

I  ! 


! 


I 
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I  I 

'  VIII.  CONCLUSION 

i  - - -  "7 

i 

The  three  methods  considered  .for  eigenvalue  problems 
yielded  satisfactory  values  for  the  eigenvalues  and  the 
eigenfunctions.  Generally,  the  relaxation  method  seemed 

J 

to.be  most  satisfactory.  It  Was  straight-forward  to  pro¬ 
gram,,  and  it  was  faster  then  Rayleigh-Ritz  and  dynamic 

I  ■  i  I 

programming.  It  converged  rapidly  even  if  step  functions 
■were  used  on  several  different  tests  figures,  such  as  the 
L-shaped  and  triangular  regions. 

1  The  dynamic  programming  method  converged  in  the  same 
number  of  iterations  as,  reliaxatlon,  but  gave  poorer  esti¬ 
mates' of  the  second  and  third  eigenvalues.  It  of  course 
was  much  more  difficult  to  program  and  required  more  com¬ 
pleting'  time  due  to  the  needed  inverses. 

The  Rayleigh-Ritz  method  seemed  to  have  little  to  re- 

i 

commend  it  due  to  the  computer  time  required.  It  required 

j 

much  finer  meshing  ip  order  to  obtain  a  satisfactory  accu- 

i  » 

racy.  The  only  advantage  it, had  was  that  no  iteration  was 
required. 

i  1  i  i 

Of  the  three  methods  considered  in  the  solution  of 
Poisson's  and  the  biharmonlc  equations  Galerkin's  method 
was  the  fastest  and  gave  accuracy  comparable  to  the  Ray¬ 
leigh-Ritz  i method .  It  was  also  the  simplest  of  the  three 
methods  used  to  program. 


i 


Dynamic  programming  gave  the  best  accuracy  generally, 
but  it  required  more  computer  time  than  Galerkin's  method. 

Rayleigh-Ritz  had  the  same  difficulties  as  it  did  in 
the  eigenvalue  problem  and  was  considered  of  little  use. 

While  dynamic  programming  required  more  time  for  the 
solution  of  both  eigenvalue  problems  and  elliptic  partial 
differential  equations,  it  was  very  powerful.  It  only 
required  a  knowledge  of  the  function  U  on  the  boundary. 

It  can  be  extended  to  irregular  regions  [l^*  It  obtained 
very  good  accuracy.  Much  of  the  time  was  spent  computing 
the  inverses.  If  the  same  points  were  used  in  solving 
several  different  problems,  these  inverses  could  be  cal¬ 
culated  once  and  thereafter  entered  as  data,  reducing  the 
computer  time  greatly. 
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COMPUTER  OUTPUT  1* ( RELAXATION) 


EIGENFUNCTIONS 

PATHs  8  OMEGAs  A. 424858 


X 

0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  *  3 

0  .  3 

0  .  3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  ,  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  ,  5 

0  .  9 

0  .  9 

0  .  9 


Y 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  ,  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 


U(X,Y) 
0.09719044 
0.2508801 
0.3094832 
0.25267961 
0.09705645 
0.2508801 
0.6512997 
0.8061690 
0.6586339 
0. 2531579 
0.3094832 
0.8061690 
1.000000 
0.8169372 
0.3137754 
0.2526796 
0.6586339 
0.8169372 
0.6658412 
0. 2551157 
0.09705645 
0.2531579 
0.3137754 
0.2551157 
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COMPUTER  OUTPUT  1 .( RELAXATION) 


61 G6NFU NCT IGN-2 

PATH=  11  OMEGAs  .6.  927361 


X 


0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  .  1 

0.3 
0.3 
0  .  3 

0  .  3 

0  .3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0.9 

0  .  9 

0.9' 


Y 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0.9 
0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0.3 
0  .  5 

0  .  7 


U  { X  » Y  ) 
0.1853049 
,  0.4886459 
0.6086538 
0.4958898 
0.1895596 
0.3032387 
0.8008109 
0. 9981067 
0.8116993 
0.3099311 
0.001278793 
0. 008393560 
0.01233941 
0  .008128799 
0.001685064 
-0.3083511 
-0.8038315 
-0.9953710 
-0.6089499 
-C. 3103257 
-0.1921543 
-0.5OC6261 
-0.619429$ 
-0.5025631 
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COMPUTER  OUTPUT  1 • ( R6L  AXATION  ) 


EIGEN FU NOT ION=3 

PATH- 11  OMEGA*  6.927361 


X 

0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  .  3 

0  .  3 

0  .  3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0  .  9 

0  .  9 

0  .  9 


Y 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  •  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 


U  ( X  *  Y ) 

0. 1857997 
0.3045443 
0.002898388 
-0. 3070407 
-0. 191655^ 
0.4894530 
0.8029368 
0.01102810 
-0. 8017028 
-0.4998162 
0.6086553 
0.9981196 
0.01235584 
-0.9953650 
-0.6194320 
0.4950783 
0.8095817 
0.005508117 
-0.  8110784 
-0.  5033802 
C. 1890559 
0.3086166 
5  .9269CG  '-0  5 
-0. 3116446 
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COMPUTER  OUTPUT  2. (RELAXATION ) 


eigenfunctions 
PATH-17  OMEGA* 


19.60767 

U(X.Y) 


0.  1011795 
0.2574610 
0.3156426 
C.  2600489 
0.1008759 
0.2580318 
0.6597930 
0.8128256 
0.6714125 
0, 2608653 
0.3149233 
0.8084638 
1.000000 
0.8271890 
0.3214785 
0.  2583916 
0.6639176 
0.8214287 
0.6778069 
0.2627081 
C. 1004035 
0.2579977 
0.3188558 
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COMPUTER  OUTPUT  2 .( RELAXATION  ) 


EIGENFUNCTIONS 

PATH= 16  OMEGA3  48.00555 


X 


Y  U(X.Y) 


0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  ♦  1 

0  •  3 

0  .  3 

0.3 
0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  7 

0  .  7 

0  .  7 

0  ,  7 

0  .  7 

0  .  9 

0  .  9 

0  .  9 


0  .  I 

0  .  3 

C  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  «  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  ,  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 


0.1909801 
0.4904021 
0. 6059132 
0.4979762 
0.1925782 
0. 3037586 
0.7884125 
0.9826701 
0.8096838 
0.3130895 
0.003429962 
0.01207077 
0.01897281 
0.01404283 
0.004551671 
-0. 3064085 
-0.7943833 
-0.9861545 
-0,8114642 
-0.3132687 
-0. 1937618 
-0.5022484 
-0.6227177 
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COMPUTER  OUTPUT  2 • ( RELAXATION ) 


EIG£NFUNCTI0N*3 

PATHs 19  OMEGA*  48.02386 


X 


Y  U(X  »Y) 


0  .  1 
0  .  1 

0  •  1 

0  .  1 

0  .  1 

0  •  3 

0  •  3 

0  .  3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0  •  9 

0  .  9 


0  .  1 
0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 


0.1851881 
0.30  20115 
0.009329475 
-0.3049613 
-0. 1935259 
0.4791670 
0. 7690626 
0.02798434 
-0.7933792 
-0.5036968 
0.5902148 
0.9785894 
0,03599443 
-0.9824019 
-0.6229793 
0.4810330 
0.7979080 
0.02484581 
-0.80  38256 
-0.5073704 
0.  1851680 
0.3065956 
0.007537059 
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COMPUTER  OUTPUT  3. (GAL ERKIN* 


VALUES  OF 

C(I*  ARE 

Cl 

s 

1.000000 

C2 

M 

0 

C3 

- 

0 

j 

X 

Y 

U(X, Y) 

0 

• 

1 

0 

• 

1 

0.008099992 

0 

• 

1 

0 

• 

3 

0.01889999 

0 

• 

1 

0 

• 

5 

0.02249999 

0 

• 

1 

0 

• 

7 

0.01890001 

0 

• 

1 

0 

• 

9 

0.008100022 

0 

• 

3 

0 

• 

1 

0.01889999 

0 

• 

3 

0 

• 

3 

0.04409999 

0 

• 

3 

0 

• 

5 

0.05249999 

0 

• 

3 

0 

• 

7 

0.04410003 

0 

• 

3 

0 

• 

9 

0.01890005 

0 

• 

5 

0 

• 

1 

0.02249999 

0 

• 

5 

0 

• 

3 

0.05249999 

0 

• 

5 

0 

• 

5 

0.06250000 

0 

• 

5 

0 

• 

7 

0.05250004 

0 

• 

5 

0 

• 

9 

0.02250007 

0 

• 

7 

0 

• 

1 

0.01890001 

0 

• 

7 

0 

• 

3 

0.C4410003 

0 

• 

7 

0 

• 

5 

0.05250004 

0 

• 

7 

0 

• 

7 

0.04410006 

0 

• 

7 

0 

• 

9 

0. 01890007 
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COMPUTER  OUTPUT  4. (RAYLEIGH) 


VALUES  OF  CU)  ARE 


Cl  * 

1.120410 

C  2  * 

-0.5808935 

C3  ® 

0.3500372 

X 

Y 

U(X,Y) 

0 

.25 

0 

.25 

0.03816108 

0 

.25 

0 

.50 

0.04937190 

0 

.25 

0 

.75 

0.03599290 

0 

.50 

C 

.25 

0.04937190 

0 

.50 

0 

.50 

0.06231650 

0 

.50 

0 

.75 

0.0446.  "59 

0 

.75 

0 

.25 

0.03599290 

0 

.75 

0 

.50 

0.04461559 

0 

.75 

0 

.75 

0.03179  573 
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COMPUTER  OUTPUT  5 . { RAYLS IGH ) 


EIGENFUNCTION*!  OMEGA* 


THE 

VALUES 

of  cm 

ARE 

Cl 

3 

29.89044 

C2 

a  0. 

003808264 

C3 

* 

2.580980 

X 

Y 

0 

.25 

0 

.25 

0 

.25 

0 

.50 

0 

.25 

0 

.75 

0 

.50 

0 

.25 

0 

.50 

0 

.50 

0 

.50 

0 

.75 

0 

.75 

0 

.25 

0 

*75 

0 

.50 

0 

.75 

0 

.75 

4*447512 


U(X,YJ 

1.073552 

1.401157 

1.028184 

1.431355 

1.868153 

1.370868 

1.073485 

1.401070 

1.028116 
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COMPUTER  OUTPUT  S.(  RAYLEIGH) 
EIGENFUNCTIONS  OMEGAs  7*230292 


THE 

VALUES 

OF  cm 

ARE 

Cl 

= 

9.664003 

C2 

= 

112.0784 

C3 

*  “112.0294 

X 

Y 

U(X,Y) 

0 

.25 

0 

.25 

0.3401793 

0 

.25 

0 

.50 

1.766417 

0 

•  25 

0 

.75 

2 .309447 

0 

o 

in 

• 

0 

.25 

-0. 8598440 

0 

.50 

0 

.  5Q 

0.  6040002 

0 

.50 

0 

.75 

1.765845 

0 

.75 

0 

.25 

-1.629948 

0 

.75 

0 

.50 

-0.8604184 

0 

.75 

0 

.75 

0.3393202 

60 


COMPUTER  OUTPUT  5  . (RAYLEIGH) 


EIGENFUNCTION-3  OMEGA-  7.247070 

THE  VALUES  OF  C( I)  ARE 


Cl 

9*656152 

C2 

** 

-112.4433 

C3 

- 

-111.6595 

X 

Y 

U(X.Y) 

0 

.25 

0 

.25 

-1.630179 

0 

.25 

0 

.50 

-0.8650630 

0 

.25 

0 

.75 

0.3325842 

0 

.50 

0 

.25 

-0.8558773 

0 

.50 

0 

.50 

0.6035087 

0 

.50 

0 

.75 

1.761141 

0 

.75 

0 

.25 

0.3463630 

0 

.75 

0 

.50 

1.770327 

0 

.75 

0 

.75 

2.309128 

COMPUTER  OUTPUT  6. (DYNAMIC) 


0.1 
0  .  1 

0  •  1 

0  •  1 

0  .  1 

0  .  3 

0*3 
0  .  3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  •  5 

0  .  5 

0  .  7 

0  •  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0  .  9 

0  .  9 

0  .  9 

0  .  9 


0 

1 

0.008099854 

0 

3 

0. 01889967 

0 

5 

0.02249958 

0 

7 

0. 01889964 

0 

9 

0.008099858 

0 

* 

1 

0.01889966 

0 

3 

0.04409914 

0 

5 

0.05249893 

0 

7 

0.04409913 

0 

9 

0.01889965 

0 

I 

0.02249958 

1 

0 

3 

0.05249896 

0 

5 

0.06249879 

0 

7 

0.05249900 

0 

9 

0.02249960 

0 

1 

C. 01889965 

0 

3 

0.04409916 

0 

5 

0.05249896 

0 

7 

0.04409920 

0 

9 

0.01889972 

0 

1 

0.  CG8099854 

0 

3 

0.01889966 

0 

5 

0.02249962 

0 

7 

0.01889972 

0 

9 

0.008099906 

i 


I 

I 

COMPUTER  OUTPUT  7.<OVNAMIC> 

I 


EIGENFUNCTION 

*1 

PATH* 

4 

1  OMEGA*, 

4.423473 

X 

i 

Y 

i  U( X. Y ) 

0 

.  1 

0  .  1 

0.09554338 

0 

.  1  1 

0  .  3 

0.2500933 

0 

.  1 

0.5 

C. 3091003 

0 

.!  1 

0  .  7 

0. 2500929 

0 

.  1 

0  .  9 

,  0.09554338 

0 

•  3 

0  .  1 

0.2500930 

1 

0 

.3i' 

(ft 

- . 

o 

0.6546414 

0  1 

•  3  ; 

0  .  5 

0.8090985 

0 

.  3 

0.7 

0.6546412 

0 

3  i  • 

!  0  .9 

0.2500929 

,  0 

.  .  5 

0  .  1 

C. 3091002 

0 

.  5 

0  .  3 

C. 8090987 

!  0 

» 

.  5 

,0  .5! 

1.000000 

0 

.  5 

0  .  7 

0.8090994 

0 

.  5 

0  i.  9 

C. 3091004 

0 

.  7 

0  .  1 

0.2500930 

0 

.  7 

0  .  3 

0.6546419 

0 

.  7 

0.5 

0. 809Q992 

0 

.  7 

0  .  7 

0.6546427 

0 

.  7 

0  .  9 

0.2500940 

0 

•'  9  i 

0  .  1 

0.09554315 

0 

.  9 

0  . ,  3 

0.2500930 

0 

9 

'0.5 

0.3091010 
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COMPUTER  OUTPUT  7. (DYNAMIC) 


EIG6NPUKCTICN 

-2 

PATH- 

6 

OMEGA- 

6*927 172 

X 

Y 

U(X»Y) 

it 

0 

• 

1 

0 

• 

1 

0.1913275 

0 

• 

1 

0 

• 

3 

0*5004815 

0 

• 

1 

0 

• 

5 

0.6183150 

0 

• 

1 

0 

• 

7 

0.5004799 

0 

• 

1 

0 

• 

9 

0.1913269 

0 

• 

3 

0 

• 

I 

0.3092604 

o 

• 

3 

0 

• 

3 

0.8089828 

0 

• 

3 

0 

• 

5 

0.9994567 

0 

3 

0 

• 

7 

0.8089805 

0 

• 

3 

0 

• 

9 

0.3092589 

0 

• 

5 

0 

• 

1 

1.959012' -06 

0 

• 

5 

0 

• 

3 

3.750642'-06 

0 

• 

5 

0 

• 

5 

2.407420*-06 

0 

• 

5 

0 

• 

7 

-4. 548319* -07 

0 

• 

5 

0 

• 

9 

-7. 790408 '-07 

0 

• 

7 

0 

• 

1 

-C.  3092566 

0 

• 

7 

0 

• 

3 

-0.8089775 

0 

• 

7 

0 

• 

5 

-0.9994535 

0 

• 

7 

0 

• 

7 

-0. 8089827 

0 

• 

7 

0 

• 

9 

-0.3092611 

0 

• 

9 

0 

• 

1 

-C. 1913257 

0 

• 

9 

0 

• 

3 

-0.5004784 

0 

• 

9 

o 

• 

5 

-0. 6183146 
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COMPUTER  OUTPUT  7. (DYNAMIC) 


E1GENFUNCTI0N=3 

PATHni  OMEGA*  6.927173 


X 


Y  U  (X  * Y  ) 


0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  •  1 

0  .  3 

0  .  3 

0  .  3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0  .  9 

0  .  9 


0  .  1 
0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 
0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 
0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 


0.1913257 
0.3092576 
-1 . 199204*-06 
-0. 3092595 
-0.1913273 
C. 5004785 
0.8089781 
-1.801975'-06 
-0.8089815 
-0.5004808 
0.6183135 
0.9994553 
1 *738089 ' -06 
-0.9994555 
-0.6183147 
0. 5004803 
0. 8089837 
4 .765  801* -0  6 
-0.8069792 
-0.5004812 
0.1913269 
0. 3092602 
2 • 62  28 19 ' -0  6 
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COMPUTER  OUTPUT  8. (DYNAMIC) 


V  U(X,Y) 


0  .  1 

0  .  1 

0  .  1 

0  .  I 

0  .  1 

C  .3 

0.3 
0  .  3 

0  ,  3 

0  .  3 

0  .  5 

0  .  5 

0  .5 

0  .  5 

0  .  5 

0  .  7 

0  ,  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0  ,  9 

0  .  9 

0  .  9 

0  •  9 


0  .  1 

0  .  3 

0  .  5 

0  .  7 

0.9 
0  •  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  ,  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 


0.008099731 
0.01889934 
0.02249917 
0.01889930 
0.008099727 
0.01889932 
0.0**409828 
0.05249788 
0.04409826 
0.01 889930 
C. 02249917 
0.05249792 
0.06249751 
0.05249795 
C. 02249918 
0. 01889931 
0.04409831 
0.05249793 
0.04409834 
0. Cl  889938 
0.008099716 
0. 01889932 
0.02249920 
0. 01889938 
0.008099772 
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COMPUTER  OUTPUT  9. (DYNAMIC) 


SIGSNFUNCTION 

=i 

PATH= 

'  3 

OMHGA- 

19,5763° 

X 

Y 

U  (  X  a  Y  ) 

0 

• 

1 

0 

• 

1 

O, 09549332 

0 

• 

1 

0 

3 

0.2500034 

0 

• 

1 

0 

* 

5 

C. 3090197 

0 

• 

1 

0 

• 

7 

0.2500030 

0 

• 

1 

0 

• 

9 

0.09549332 

0 

« 

3 

0 

• 

1 

0.2500032 

0 

• 

3 

0 

• 

3 

0.6545131 

0 

3 

0 

• 

5 

0.8090189 

0 

• 

3 

0 

• 

7 

0.6545127 

0 

• 

3 

0 

• 

9 

0.2500032 

0 

• 

5 

0 

• 

l 

0.3090197 

0 

• 

5 

0 

• 

3 

0.8090194 

0 

• 

5 

0 

• 

5 

l. 000000 

0 

« 

5 

0 

• 

7 

C. 8090203 

0 

• 

5 

0 

• 

9 

0.3090200 

0 

• 

7 

0 

• 

1 

0.2500033 

0 

* 

7 

0 

• 

3 

0.6545139 

0 

• 

7 

0 

* 

5 

0.8090  203 

0 

• 

7 

0 

• 

7 

0. 6545147 

0 

• 

7 

0 

• 

9 

0.2500044 

0 

• 

9 

0 

• 

I 

0.09549326 

0 

♦ 

9 

0 

# 

3 

0.2500034 

0 

• 

9 

0 

• 

5 

0.3090  20  6 
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COMPUTER  OUTPUT  9. (DYNAMIC) 


EIGENFUNCTIONS 

PATH-  4  OMEGA-  47.98405 


X 


U(X,Y) 


0  .  1 

0  ,  l 

0  .  1 

0  ,  1 

0  .  1 

0  .  3 

0.3 
0  .  3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .5 

0  .  5 

0  .  5 

0  o  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  9 

0  .  9 

0  .  9 


0 

« 

1 

C.  1911G72 

0 

3 

0.5001713 

0 

5 

0.6181207 

0 

7 

C.  5001690 

0 

9 

0.1911063 

0 

1 

0.3091234 

0 

3 

0.8090394 

0 

5 

0.9998273 

0 

7 

0.8090362 

0 

9 

0.3091207 

0 

1 

“9, 384042 1 “07 

0 

3 

8,968278* -06 

0 

5 

l .008165*-05 

0 

7 

3. 054505* -06 

0 

9 

“4  ,620023'-06 

0 

1 

“C. 3091281 

0 

3 

-0.8090287 

0 

5 

-0.9998166 

0 

7 

-C.  8090367 

c 

9 

-0. 3091323 

0 

1 

-0. 1911175 

0 

3 

-0.5001757 

0 

5 

-0.6181255 
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COMPUTER  OUTPUT  9* (DYNAMIC) 


El GEN FUNCTION* 3 
PATH-  7  OMEGA* 


47.98486 

U(X*Y ) 

-/ 

0.1911126 
0.3091245 
-9 .687 110 *“07 
-0.3091258 
-0.1911120 
C.  5001729 
0.8090308 
-2.217028*-06 
-0.8090366 
-0.5001741 
0.6181221 
0.9998204 
1.343255*-06 
-0.9998225 
-0.6181232 
C. 5001751 
0.8090382 
5.371387'-06 
-0.8090329 
-0.5001737 
0.1911140 
0.3091279 
3.674680'-06 


COMPUTER  OUTPUT  1 0. ( GALERKIN) 


VALUES  OP  C(I)  ARE 

Cl  «  1.596382 
C2  =  1.596392 
C3  =  -2.598899 


X 


Y  UtX.Y) 


0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  .  1 

0  •  3 

0  .  3 

0  .3 

0  .  3 

0  .  3 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  5 

0  .  7 

0  .  7 

0  .  7 

0  .  7 

0  .  7 


0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 

0  .  1 

0  .  3 

0  .  5 

0  .  7 

0  .  9 


0.002375632 
0.01059511 
0.01862749 
0.02069907 
0.01103618 
0.01059507 
0.03192532 
0.04658196 
0.04633235 
0.02294400 
0.01862741 
0.04658184 
0.05916639 
0.05281764 
0.02397245 
0.02069896 
0. 04633217 
0.05281758 
0.04240138 
0. 01732974 
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COMPUTER  PROGRAM  1. 


RELAXATION  METHOD  APPLIED  TO 
VaU  •  -X2U,  EQ.  (1.1) 


BrGIN 

P™OSFIGi  nvAUF  S  AND  cl GCf.'V?  ctors  by 
P=L|XAT!0N^BAS2D  ON  RAYLEIGH'S  FORMULA  ON  ANY 

FOP  MALPARAMei?  RS: 

OM  OM’GA 

CM2  DM  IGA  SQUARED 

m:  ksm.®  of  nM_GA  soimpeo 

N1  M  AX.  NUMBER  OF  X  POSITIONS 

°1  MAX.  NUMBER  rc  Y  POSITIONS 

HI  MESH  OF  X  VA^IiBLP 

H2  MESH  OF  Y  VA.eiABLF 

N2  ARRAY  OF  LOW-R  POSITIONS  OF  Y 

N3  ARRAY  OF  U0!:,rR  DOSITIONS  OF  Y 

U,U1,Z  UNKNOWN  E  I  G  'NF  UNCTIONS 

P.F  POTENTIAL  "NT?.  GY 

Kl  KINETIC  5iN’:P«Y{ 

RPAL  K I  •  P': . OM 2 tOMO 2 «  ,  Z  StOM.ZNORM; 

RPAL  H.H3.H4;’ 

OM2s=5COO.^; 

Ns=  l; 

~Rt :=C .002 ; 

CAT  1:  =  50C0, 0 ; 

H 3  S  =H  1  •*  <2; 

H4:  =W2  ;  '"2 ; 

*'r  r- 1  N 

R?AL  ARRAY  X(0::N1); 

RIAL  ARRAY  Y(0;:?l)i 
INTG'S  ARRAY  N2  »N3  (  * : : Ml  ) ; 

RF AL  ARRAY  Z.lU'JKc:  :N1.0::P1)  ; 

FOR  I :  =''■  U  NT  T  L  N1  ro  * 

FOP  J : =0  UNTT  L  PI  DO 
BEGIN  UK  I .  J  | ;  =c  . 0  ; 

Z( I j  J )  :»c .0; 

END ;  ? 

FOR  !:  =0  UMT'L  mi  do  X ( I  ) :  =  t  -H  1 ; 

pop  !  ;=n  UNTIL  PI  DO  Y(I):  =  I>H2? 

FOR  I  :  =0  IJNT?L  NI  OP  R  "  APPU  (  N2  ( I  )  )  J 

fqf  I:=0  UNTIL  NI  DO  ?.  ;  A  DON  ( N3  (  I  ))  ; 

IF  L=1  THVN 
BEGIN  1:=?; 

=0R  T  :=r  ‘JNT  T  l  NI  O'” 

F0°  J  :  =  f.  UN  Til  5 1  r-n 

r!o!fr,i:(X(I,"X(n'  2  1  <Y(J)-Y(J  )*■  2); 

;M0; 

IF  L=2  Tw^M 
R*  GIN 

S0R  I  :=  3  UNTIL  Ml  OP 
FQR  J :  =r  UNTIL  =>1  "0 
BOGIN 

U(I»J):=Z(I»J  >; 

Z(I  t  J  )  :  =  CX  { I  )-:<  (I)  2  )"(Y(J)-Y(jrU2)MC 

„  _  END; 

Cat :=K1 ;L :=3;N:=1;0M2: =5000.0; 


(  Y  (  J  )  -  Y  {  J  )* 


(Y(J  )  -Y ( J  )  *  *  2  M (0. 5- X  (  I)  ) 


71 


GO  TO  Di 
“ND ; 

IF  1*3  TH'N 
BF6IN 

f6R  1**0  UNTIL  N1  DO 
FOR  JS *0  UNTIL  PI  DO 
BEGIN 

Ultl » J) :»Z!I T  J>  ? 

ZU»J):  =  (Xm-Xl!  )  "*‘*2 )  *r  { Y  ( J  )-Y<  J)  **2)  (0.5-Y(  J)  )  ? 
END? 

CATl:*Kf~  ?L:=4?N:=l?OM2:  =  5000.0t 
GO  TO  D? 

END  ? 

COMMENT  GET  ?=  AND  <c? 

A:P2s*K?s=0#0? 

FOR  I s  =0  UNTTL  Nl-1  DO 
FOR  J:*N2(I)  UNTIL  N3<n-1  DO 
o“GIN 

PE :*P~+ !(!Z<I+1»J)-ZC I '  2  )  H 
♦  ((lZ(I,J+ll-2n.JI  )  /H  2  )  •’  Z)  ■  H; 

K • s  * KE+ ( Z ( ! , J I *  2 )  H 1 * H2  J 
'NO; 

0MG2:=0M2; 

OM2:=Pr-/Kf:?OMJ=SORT(OM2)  J 
TF  (0M2>0M02)  OR  (N>40)  THEN 
3  “GIN 

WRIT” ("NOT  CONVERGING” )? 

WRI  TE  <  ”M =" ) ;  w? !  TE  CM  ( N ) ; 

WRITE! ”0^02  AND  3^2  5V»FIT?0N(0M02*0M?»  ; 

GO  TO  R; 

END : 

IP  APS (0M2-0MD2  KTRA  THEN 
BEGIN 

IOCOUTROL! 3) ? 

W  R  T  T  ”);WRIT-(”  ")?WPTTM" 

WFITt"  COMPUTE P  OUTPUT  1 .( RELAXATION  ) 

")  ? 

WRITE!”  «);WRITC(”  » ) ; 

I NTFI ELDSI  Zv  :  =1 ; 

WRIT,-!”  -IGHNrUNCTI0N=".L-U  ;WRIT”(” 

I N T  c  !  "  L  DS I  Z  E  s  =2  ? 

WRITE!”  FATH  =  ”,N,”  OM~GA=”,OM); 

WRIT;:!"  ”)  ;WP  ITr!”  ”  )  ? 

WRIT;!”  X  Y  ”, 

”  U  (  X  » Y )  ” )  ? 

WRIT*!”  ”); 

C0F.  I  :  =  1  ST'*  ?.  UNTIL  M  -i  DO 
FOR  J:  =  1  STE?  2  UNTTi_  pi-i  on 
BEGIN  !  MT  c  I  EL  DS I  Z  - :  =2  ;  " 

WRIT”!”  ”,I  DIV  N1 REM  Nl,«  ”, 

j  njv  ®  1  »***.*•  J  R“M  Pl.Z(I.J)); 

WRIT?!"  ” ) ; 

CND  ?~ND  ? 

IF  C  A5S  (0M2-nMr2  )  <"  R  A  )  AMO  !L<4)  TW'N  C-0  TO  e? 

IF  -*,BS(0M2-0MJ2)<;  RA  THEN  GO  TO  R; 

N  :  =  N+ 1 ; 

COMMENT  GET  ”"W  Z; 
for  I : =1  UNTTL  Nl-1  00 
BEGIN 

=0P  J  :  =N2  (  I )  +1  UNTIL  NSUI-l  DO 

Z<1  :  =  <N4v(Z(I  +  l,J)+ZU-i  .  J  )  )  +H2  ,;  (  Z  ( I  ,  J+ 1  )  +  Z  ( I ,  J  - 1 )  )  ) 

/!2.‘'H3+2.  H4-^M 2  ‘ (  H •i-  :'2)  )  ? 
rND; 

IF  L>2  TH'N 
BEGIN 

Cl :=DOG:=OOG1 :*C.C ? 

*=0R  I  s  =1  UNTIL  Nl-1  CO 
=0P  J :  =  M2 (  I )  +  1  UNTIL  M 3 ( I  1-1  CO 
BOGIN 

C 1 :  =C  1  +  Z ! !  » J  J  m; 

OOG:=00G+Zi I, J )  U( I , J )  Hs 
DCGL :=OOG1  +  Z( I , J)  U1  (  I ,  J  J  H ; 
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MP  ? 

C2 : =OOG/P  ^  ;C3: ^POOl/C/Tl ; 

GriC'  I  :=1  !‘MT*L  M-l  0^ 

FCR  J:  sf.‘2(  T  )  +  l  UNTIL  N3(!)-l  00 
.  Z  ( T  *  J )  :  =  Z(I»J  )“C1~C2V‘U(  I  *  J  I  “C  3*U  1  ( I 
MP  5 

n  •  .i  • 

f=r»  U-(,  UNTIL  Ml  DO 
C0R  J  :  =  N2  ( !  )  UNTIL  N3< !  )  DP 
T-GI‘1 

zss«abs(Z(  i».i ) ); 

Ic  ZS>ZMO?M  TH  N  ZNORMjsZS; 

r  Mr; 

cos  I: -A  Ijh.7 

rC’F  J:=NI2 <!) 
on  to  i,; 
if  MO; 

R : :  NO. 


'  L  M1  r'"' 
i  IMTI  L  N3  ( I  ) 


00  Z(I ,J) : =Z ( I » 


,  J  »  S 


J ) / ZNQRM 
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COMPUTER  PROGRAM  2. 
RELAXATION  METHOD  APPLIED  TO 
*  XaU»  EQ.  (1.3) 


ZNQRM 

H 

UfOl.Z 
PS 
KS 


itgin 

CCWVsMT  eINDS  sIG'NVALUSS  AND  SIG-NV“CTOPS  BY 
R'LAXATTOM  BASED  ON  PAYL'IGH'S  FORMULA  FDR 
0  4 { U  )=W '  *2  U. 

FORMAL  PAR  AMS  T'  PS  : 

CM  OM*fiA 

CM2  OMEGA  SQUAR'D 

OM02  LAST  VMIJr  D=  QM“GA  S  QUARTO 

ZNOPM  MAX.  NOPM 

MESH  SIZE  OF  X  AND  Y 
UNKNOWN  £  I G7NC UNCTIONS 
POTENTIAL  NS  RGY 
KINETIC  EM’rPGY; 

REAL  ARRAY  Z,U,U1 (C: : 14,0s : 14  )  S  „  , 

REAL  Kg.PL:.OM2.0M02.:;RA  .  ZS.OM.ZNORM.H.CAT.DOG.CATl  .DOGli 
REAL  Cl,C2»C3* 

INTEGER  N’t L ? 

REAL  ARRAY  Y(0::10>; 

REAL  ARRAY  X (0 « slO  » S 
LOGICAL  SWITCH; 

REAL  OMO; 

0^1=7000.0*, 

SWITCHs=TRUP? 

CAT  :-CAT I : =2000.0 ; 

H ; =0.  l  •  °M2 : --2000.0;  n  :  =1 ;  ?R  A:  -0 .003 ;  L :  =  1 ; 

OM 2  *.  =  700000.0; 

ERA:  =0  .005; 

FOR  I  :  =0  UNTIL  10  DO  Xt  I)  :  =  YU  )  J  =  I*H; 

CCR  I : =0  UNTIL  14  OD 

FOR  J;=C  UNTIL  14  DO  Z  U  ,  J  ) :  =U  U  ,  J  ) :  =11 1 U ,  J  ) :  =0 . 0  { 

B'-G  IN 

PROCEDURE  NORMA* 

3 EG IN  Z NORM: =0.0; 

FOR  I: =3  UNTIL  11  00 
FCR  J : =3  UNTIL  11  DD 

BEGIN  ZS : =ABS (Z(I *  J  M  ; 

If  ZS>ZNORM  THh  N  ZN  3RM :  =  z  $ ;  END; 

FOR  I :  =0  UNTIL  14  00 

FOR  J : =0  UNTIL  14  DO  Z ( I  *  J  ) : =2 (I , J ) / ZNORM; 

GO  TO  A; 

END  MOR'-'A; 

PPCCEDURO  FiJNCT; 

BEGIN 

IF  L  =1  THEN 
BEGIN  L : =2; 

FOR  I : =0  UNTIL  10  DO 
FOR  J :  =0  UNTIL  10  00 

Z(I+2*J  +  2):=(Xn)-XiI)**2)MY(J)-Y(J)*T2); 

NORMA ;2ND; 

IF  L  =  2  THr- Kl 

BEGIN  L  :  =  3  ;C  AT jNs  =  1*,0M2:  =  2OOO.  0; 

CM: =7000.0;  ,  , 

FOR  I;-0  UNTIL  14  DO  FOR  J:=0  UNTIL  14  DO  U  (I « J  ) : - Z ( I , J ) ; 
FOR  I : » 0  UNTIL  14  DO  FOR  J : =0  UNTIL  14  DO  Z(I.J):=0.0; 

FOR  2  s -0  UNTIL  ID  DO 
FOR  j:=0  UNTIL  10  DO 

Z(I+2.J  +  2»S  =  (X(I )-X( I )**2)*( Y( J )-Y( .5-X(I  n ; 
NORMA; END; 

IF  L=3  THEN 

BEGIN  L :  =  4 ; C AT  1 :  =  <  T: ; N :  =  1  *, 0 M 2 :  =  2 0 0 0 . 0 *, 

OM; =7000.0; 

FOR  I  :  =  0  UNTIL  14  DO  C0P.  J:=0  UNTIL  14  DO  U1  U . Jl :=Z ( I » J ) 
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B0R  1:@0  UNTIL  14  00  FOR  J:=0  UNTIL  14  DO  Z*I,J):=0.0{ 
PHR  I  :=0  UNTIL  10  00 
FOR  J:*=0  UNTIL  10  DO 

Z  ( 1+2  t  J+2  ):  -  (XC I  l-Xim*2WylJ)-YIJI^2IM.5»Y(J|]  ; 

normajsnd; 

END  FUNCT; 

IF  SWITCH  THEN  BEGIN  SWITCH: -FALS S {FUNCT ;2NDi 
COMMENT  GST  PS  AND  KS { 

A:P5:=K2:=0.0; 

FOR  1**2  UNTIL  1 2  DO 
FOR  J:  =2  UNTIL  12  DO 

BEGIN  PS: -PS+ ( Z  < I  +  1 , J )-4. *Z( 1 , J 1 +Z ( 1- 1 , J ) +Z (I , J+ 1 ) 

+Z<  I  , J-l) ]**2/*H*H){ 

KS:=K£+(Z<  I.J  )=vH)**2;BND; 

OMO :  -OM; 

OH2*  =  PE/KS;OM:aSORT(OM2) ; 

IF  (OM>OMO )  OR  <N>60)  THIN 

eSGIN  WRITE* "NOT  CONVERGING  ">{ 

WRITS*"  ")  {WRITE* »N,ZNQRM,OMO,OM») { 

WRITE (N» ZNORM ,0M0 ♦ OM ) {WRITE ("  " ) { 

FOR  I : -1  UNTT  L  11  DO 

WRITE* Z( 3i I ) ,Z< 5.  I ) • Z ( 7. 1 > .Z*9.I> .Z(11.I >){ 

GO  TO  R ‘END; 

IF  ABS*CMO-OM)<FRA  THEN 
BEGIN  IOCOMTROL (3 ) { 

WRITE*"  WRITE*"  "I {WRITE*"  ">{ 

WRITE*"  COMPUTER  OUTPUT  2.", 

" (RELAXATION) "> { 

WRITE*"  "){ WRITS* "  ") { 

INTFISLDSIZE : =1* 

WRITS*"  5IGHNFUNCTI0N=«tL-l)  {WRITE*  «  '•)? 

INTFIfLDSI ZS : =2 { 

WRITE*"  PATH=".N."  OMSGA=" . GM  I ; 

WRITE*"  ") {WRITE*"  ") { 

WRITE*"  X  Y  ", 

«  U  *  X  ♦  Y ) "  )  { WP.  I  Tc  (  "  ")  {WRITE*"  »){ 

FOR  I :  =1  STEP  2  U Nr  I L  <3  DO 
FOR  j:-l  STSP  2  UNTIL  9  DO 
BEGIN  INTFI SL  OSIZ  S  s  =2 • 

WRITE*"  «,I  DIV  10,".", I  REM  10, 

"  «,J  OIV  10,".", J  ROM  10  »  Z(  1+2  ,  J+2  )  )  { 

WRITS*"  "){ 

END; 

IF  L<4  THEN  FUNCT  ELSE  GO  TO  R{ 

END; 

N : *N+1 5 

COMMENT  GET  NSW  Z{ 

FOR  I: =3  UNTIL  11  DO 
C0R  J: =3  UNTIL  11  DO 

2  ( I  *  J  I  :  =  *  R.vZ*  1+1  ,J)t8.*Z*I-l,JI  +  8,*Z<l,J  +  l) 

+8.---Z*  I  ,  J-l  )-Z*!+2,J>-Z(I,J-2)-2.*Z*I+l,J+l) 

-2.  +  Z(I-l  ,J-1 I -2 .  >  Z  < I -1 ,J+1 )-2.*Z*I+l,J-l ) 

-Z(  1-2,  J  >-Z  (I .  J  +  2  1  )  /  (  20.  -0M2*  Hr  H*  H*  H )  { 

FOR  I:  -2  UNTIL  12  DO 

BEGIN  2*0,1) :=-Z(^, I) ; Z(14, I) :=-Z(10,I) { 

Z(l. I  ):=-Z<  3. I ) ;Z  < 13. 1 ) :=-Z< 11.! > { 

Z(I  ,OI:s-Z(!,4) ;Z( I ,14):  =  -z* I, 1C ){ 

Z ( I  ,1) :=-Z< 1,3) ;Z(I ,13) :=-Z*I ,11) {END; 

IF  L>2  THEM 

BEGIN  Cl :=OOG:=OOG1  :=0.0{ 

FOR  1**3  UNTIL  11  DO 
FOR  J :  =3  UNTIL  11  DO 

BEGIN  Cl:=Cl+Z(I,J)*H+H; 

DOG:  sCnG+Z  (  I,  J  )  “=  U  *  I 
DOGl:=DOGl+Z*  I,  J  )*U!  (I,  J  )»<H*H;ENDS 
C2: sDOG/CA  T; C  3 :=D0G1 /CAT1 ; 

FOR  1 :  =  3  UNTIL  11  DO 
FOR  J : “3  UNTIL  11  00 

Z  ( I , J ) :  =  Z (  I  , J»-C 1-C2"U*I , J> “C3+U1 *  I  ,  J  )  { 

FOR  I :=2  UNT I L  12  DO 
BEGIN 

Z(0,I )  :=-Z*4,n  ;z*  14. 1)  :=-Z(10.I  >;z*l  .1  )  :"-Z(3  ,I> 
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1N0{ 

IF  L<5  THEN  NOP.MA;  ,  • 

5ND ;  1 

p  ;rtjD. 


s-zn  .io) 
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COMPUTER  PROGRAM  3. 

*  GALE&KIN ' S  METHOD  APPLIED  TO 
7*U  »■  2(x2+y2“X-y ) „  EQ.  (1.5) 


COMMENT  SOLVES  PARTIAL  .DIFFERENTIAL  EQUATIONS  OP 
THE  FORM  L  ( U, )  =  F  BY  GAl!sRKIN‘S  METHOD  ON  ANY  DOMAIN. 
FORMAL  PARAMETERS i 
N  NUMBER  OF  EQUATIONS 

•  M2  MAX.  NUMBER  OF  X  POSITIONS 

M3  MAX.  NUMBER  OF  Y  POSITIONS 

N2  ARRAY' OF  LOWER  Y  POSITIONS 

K3  ARRAY  Oc  UPPER  Y  POSITIONS 

F  KNOWN  FUNCTION  . 

U  UNKNOWN  FUNCTION 

;  C < I  )  ESTIMATING  FUNCTIONS 
Dim.  L<  2STI  MATING  FUNCTIONS) 

Ad)  CO*F.  OF  ESTIMATING  FUNCTIONS; 

REAL  H  . S ,DOG , CAT , OE V , AMUL , T , OOM ; 

REAL  HI; 

ROA-L  H2.H3; 

INTEGER  M2,M3; 

INTEGER  N,Mt,R; 

E;  =0.000001  *,P:  =0  ; DEV  :  =1.0; 

READ  ( N  » M2  ♦  M3  ,H2*H3.)s 
1  Hls=H2’KH3; 

BEGIN  1 

REAL  ARRAY  X(0::M2); 

REAL  ARRAY  Y<  02  s M 3  )'; 

REAL  ARRAY  fi  (  1 :  :  N ,  1 :  *.  N  +1  )  ; 

REAL  ARRAY  A  * C  C 1 2  :N) ; 

REAL  ARRAY  r ♦ U ( 0 i : M2 . 0 : : M3) ; 

REAL  ARRAY  D 1 01 ( 1 : :N tO  : : M2 ,0 : : M3 ) { 

INTEGER  AR»AY  N2 ,N3< 0 : :M2 ) ; 

FOR  l:=0  UNTIL'  M2  00  X(I):  =  I*H2; 

FOR  T*=0  UNTIL  M3  DO  Y!I):=I*H2; 

FOP  18*0  UNTIL  M2  CO  REAOON! N2!I )) ; 

FOR  I  a  *0  'UNTIL  M2  DO  READQN(N3<  I ))  ; 

FOP'  I  :  =0  UNTI  L  M2  DO 
FOR  JS=0  UNTIL  M 3  DO  ' 

BEGIN 

f ( f  * ji  2  =0. 0 ; 

U!I  .J ) :  =  0 .0  ? 

FOR  L:=l  UNTIL  N  00 

C  (L » I  » J  )  :  =01  (  Lt  I  ♦  J)  ;  =0. 0* 

B?GfN  Comment  solve  here;1 

PROCEDURE  WRITER; 

BEGIN  WRIT  (  "S  IN  GULAK  M  )  ;  GO  TO  0; 

END  WRII^a; 

PROCEDURE  PERFORM; 

BEGIN  COMMENT  GET  U!UJ>  AND  A(I>; 

WRITE!  "  |  WRITS!"  ")*,  WRITE!"  "); 

WRITE! "  - 

WRITE!"  ")  {WRITS ("  " ) ; 

WRITE!" 

WRITE!"  ",); 

FOR  I ; =1  UNTIL  N  DO 
BEGIN 

All  ) :  =  B( I . N+l ) ; 

INTFIcLDSl Z  E  2  =1  ; 
wo  ITS!" 

END ; 

FOR  J!=0  UNTIL  M3 -1  CO 
BEGIN 

FOR  I : =0  UNTIL  M2-1  CO 
BEGIN  DOG;=O.Ci' 


COMPUTER  OUTPUT  3 . ( GALERK I N) " ) 
VALUES  OF  C ( I )  ARE" )» WRITS ( "  "); 


C'M  » "=" t A ( I ) ); WRITS!"  «); 
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FOR  L:  =  l  UNTIL  N  00  OOG t=DOG+A( LI*D(L.I .  J) ; 

U(I  ,J)i=DOG; 

S NO i END  5 

WRITER"  *  )  JWR I TS  (  "  "  WR I TE (  "  ")  ; 

WRITE <"  X  V  ». 

"UtX.Yl  '•)  5 

WRITE ");WRITS<"  " ); 

FOR  l:«l  STEP  2  UNTIL  M2-1  DO 
FOR  J:  =  1  STEP  2  UNTIL  M3-1  DO 
BEGIN  INTFISL0SI2 E: =2 { 

WRITE!"  «,I  DIV  M2. ".".I  P.EM  M2» 

"  «tJ  DIV  M3  •♦*«•'*  J  REM  M3.UU.jn; 

WRITE!"  » » ; 

END: 

GO  TO  u; 

END  PERcORM; 

PROCEDURE  SWITCH; 

BEGIN  COMMENT  SOLVES  FOR  A ( 1 )  COcF F  ; 

INTEGER  N1  ,  Ml  ; 

Ml :  =M; N1 : =N+ 1 ; 

FOR  I*«l  UNTIL  Ml  DO  B(  I  .Ml  +  1 » :*C  (  II  j 
FOP  K; =1  UNTIL  Ml  DO 

BEGIN  IF  K=M1  THEN  GO  TO  P; 

FOR  I:=K+1  UNTIL  Ml  DO 

BEGIN  IF  ABS(B!K,K))<ABS(B(I.K))  THEN 
BEGIN  R:»R+l: 

FOR  J;=l  UNTIL  N1  DO 

BEGIN  T:=S(I.J);B!I»J);=8!K,J);B!K»J):=T ;«ND 

END  ;ENO ; 

P:IF  ABS  (  8  (  K.K  IKE  THEN  WRITER  ELSE 

BEGI N  DE  V :  =DH  V+  8  <  K  ,  K  ) ;  OOM  :  =  B  <  K  ,  K  > ; 

FOR  J : - 1  UNTIL  N1  DO  B!«,J) :=B(K.JJ/DOM; 

FOR  1 i =1  UNTIL  Ml  CO 

BEGIN  AMUL:=B( I ,K) ; 

IF  I  =  K  THEN  ELSE 
BEGIN 

FOR  J :  =1  UNTIL  N1  DO 

BEGIN  8 ( I  .  J  )  :  s8  (  I  » J ) -AMUL*B ( K . J)  ; 

END; END; END?  END 5 END; 

PERCQFM; 

END  SWITCH; 

COMMENT  GUTSSSD  RJNCTIONS  D ! I  I  HERE; 

FOR  J :  =0  UNTIL  M2  00 
FOR  K:=N2(J)  UNTIL  N3(J)  DO 
BEGIN 

F(J«K):=2.*(X(J)**2+Y(K)**2-Y(K)-X(J)); 

D(  l.J.Y  ) :  =  (  X(  J  )-X<J  )*  *2  )■*<  Y(  K  )  —  Y  <  K  )«  +  2)  ; 

C(2»J»K):  =  (X(J)**2-X(J)S,*:3)1MY!K  2-Y(K  )Y=i>3  >  ; 

Cl!  It  J  »K)  :  =  2.*  <  X!  J)**2+Y(  KH*2-X(  J)  -Y  (K>); 

CH2.J.K):  =  (2.-6,^(JI)«(Y(K  )*■*  2-Y(K  )**  3) 

U2.-6.n(KH'  (X  {  J|**2-X(J 

END; 

COMMENT  SOLVE  FOR  INTEGRALS  F*0  AND  L  D*D » 

FOR  I  S  =1  UNTIL  N  00 
BEGIN  D0G:«0.0; 

FOR  j; =1  UNTIL  M2  DO 

FOR  K: =N2( J) + 1  UNTIL  NB ( J )  DO 

COG:=DOG+( F(J ,K l*0(I,J.KI*Hll; 

C ! I ) :*DOGS 
FOR  J :  =1  UNTIL  N  DO 
BEGIN  CAT :  =  0. C ; 

FOR  K :  =1  UNTIL  M2  DO 

FOR  L:-N2!K)+1  UNTIL  N3(K)  DO 

CAT :=CAT+(C1(J«K,L)AD( I.K.L )*H1) ; 

8(1 . J ) : =C AT; 

END; 

IF  I=N  THEN  SWITCH; 

END; 

END ; END; 

Q: END. 
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COMPUTER  PROGRAM  1| . 
RAYLEIGH-RIT2  METHOD  APPLIED  TO 
V2U  ■  2 (x2+y 2-x-y ) ,  EQ.  (1.1) 


P--GTN 

V:\iT  $t(.  \!  S 

P  !i  YL  1  GH-R  IT  z 


0  U  T!  f  L 
M  TM30. 


0!B|!?R:KmL  FOUATI QNS  BY 


Cl?  V  * 

L  *fiz- 

M  -  T  ’  ^S* 

H2 

OX 

*r.  t  NO 

WA 

OY 

■■  C  IN*R 

M 

numd 

P  .-ic  OOUATIONS 

C 

KNOW 

M 

FUH  TI 

j  M 

m? 

URP  ' 

LIMIT 

OF  X 

V3 

IIP  P 

R 

LIMIT 

OF  Y 

U 

UN  WN 

WN  rUNf 

71  ON 

NO 

URP" 

A 

BOUNDS 

n  m  y 

N2 

tow 

c 

BOU NOS 

ON  Y 

A(  I  ) 

CO'  = 

0=  ?STI MATING 

Mil 

-ST! 

M 

(IT! N* G  PUNOT  IONS 

F UNC  T I  ON  S 


■  t  < 


R*AL  H.i, 

INT "Gr°  M i.  *  - 
-jsC.rooooi  :F. t  « 5 ; r ;y 
R  T A  D  U  i » '•I  2 »  M3  ,  H  2 » H  3 ) » 
Hli =M2-W3; 

BOGIN 

(UAL  ARRAY  X  (  C  s  t  v'  2  > ; 
RFAL  t  FRAY  Y<o:jM3); 
PivAL  ARRAY  R  ( 1 : :  N ,  1 : 
RZ-AL  ARRAY  A  *  C  ( 1 :  t  N  ) 

R  '  AL  ARR/v  c  . - 

R'-AL  ARRAY 
INTcGTR  A*1 


:M  +  1  ) 


,U <0  ssM2 


0,01,02(1  t  :U,0::M2fC:  i  M2  »  ; 
AY  N!2  ,U3(  C  i  :M2  ) 


YU  )  :*1  <F3; 
k;  AT?M(N2<  I)  )  ! 

M2  (ill; 


02  (  l  » 1  ,  J  ) :  =0  «0 


rJ  ) 


FCR  1*«0  UNTIL  M2  Dn 
POP  1:=A  i lN!T!  L  M3  DO 
P0D  I : =0  1 1  NT i L  M2  00 
=oq  1:=-'  ijNTTt  m'  nn 
POP  I  i  =0  UNTIL  M2  DO 
FOR  J :  =0  UNTIL  M3  DO 
Pp  GIN 

FC I ,J) : =U< I  , J  > i =0. 0 
=nc  l : = 1  UNTIL  N  "O 
0  (  L  , : , >))  :  *0  1  (  L  , ! ,  j ) 

-*■10  * 

FOP  It  =  1  UNT  !  L  M  or> 

FCR  J:  =1  UNTIL  N  +  l  09  6(1 
B?GTN 
COMMrNT  SOI  V.  H 'T R  ; 

PROCHPiJR'I  writ.  Ri 

9?G!M  WR I T " ( »S IMGULAP" ) ;G0  TO  Q; 

r  NO  WPIT-R; 

PPOC'DIJP"  F‘.r<cORM; 

9"  G  IN 

CGMM;mt  0  :T  IJU  ,J  > 

WR I  T  -"  (  ••  <n  •  u?  T  T’..  I  n 
W  R I  T  U  >• 

WR !  r  (  « 
wp  1  Tv  <  '• 

C0R  I : =1  UNTIL  H  CO 
BEG  I  N 

A  ( !  )  :=B( 1 ,M+1) ; 

!MT  FIELDS  IZO:  =  1  *, 

WR!  T  ('» 

TND ; 

FOR  I:  =1  U'!T  ’l  **?.-! 
cr R  j:=i  untti_  M3-1 
Q  - GIN  DnG : =0. 0 ; 

FOR  L: =1  UNTIL  N 


0  :  T  IJU„I 
")  ?  WP  I  TP  ( 

") ;  W~ I T-  ( 

UNTIL  N 


AM  C  A  <  I  ) 
") ; WRI T 

" ) ; wr i r  ; 


< "  " ) ; 

CO^PUT-P 

< "  "  ) 


OUTPUT  k*  (  °A  YLP  I  GH  )  »  ) 


VA  L  U :  S  OF  C(I)  ARi"  )  ?  WP.  I  T  '  <  "  "); 


C",I,"=",A(I)l;WP.IT'(«  '•) 


D 

r  1 


00  DOR:  =  DOG+A ( L }"D ( L , I , J ) 


U(!  VJ) i=0DG; 

•ND ; 

WRJT|j«  "UWRITrt”  w  )  J W 3  1  T  ■  ( '•  <•); 

U(x,Y)«);  X  Y  "* 

WRITE!"  ")  jWPITf-C  ••  ); 

POP.  1:=^  5T"P  10  l  IN  TT  L  f'2  DO 
P°R,J£t0  ST”P  1*3  UNTIL  M3  00 
BEGIN  I  NT  c  I  ?L  DS  I Z .  =  2  ; 

WRI  T...  <  '•  «,I  oj»/  M2  ♦  “.“f  1*100  PgM  M2, 

URJj,(„  -!V  M-» "•"» J" loo  f.;:m  M3;6(i,j))i  41 

"NC  I 

GO  TO  0; 

END  P? PFnRM ; 

PROCHLUP'-  SWITCH? 

BEGIN 

COMMENT  SOL  V  S  FOP  A  { T  )  f;CcF'=J 
I  NT  “G'  P.  Nl.Ml; 

MU  =N;  N1  :  =  M  +  l  ; 

FOP.  1  :  =1  Until  Ml  DO  ?  (I  ,  mi  +  i  )  ;=C  ( I  >  ; 

FOP  k:=i  UNTIL  Ml  DO 
BYGIN 

IF  K  -M  1  TH “ N  GC  TO  P; 

FOP  I  J=K+1  UNTIL  Ml  00 
BOGIN 

IF  ABS ( R ( K, K ) )<ABS(B(I«K) )  THIN 
BEGIN 

R  t  sp  4- 1  • 

FOR  J:=l  UNTIL  Ml  00 

=NO;*-nd- 3?r,IN  Ts  =  en*J>  *B(I,J)j«B<K,J1  ;B  <  K , .J ) :  =T  ;^D 

P!  IF  ABS!  °  <  K  ; K  )  )<:  TH'-. N  UP  I  T„  3  r-LSii 
B  ”  G 1 N 

D~  V :  *00  W 13  (  k  ,  K  )  ;  0  DM  s  =  g  ( i/ 1  k  ) ; 

Trl  i\U  H3tI[  M  K  *  J )  J  =B  (  K  ♦  J )  /DOM; 

B"0!N 

A “  1  Jl  :  =3 (I » K  ) ; 

1“  I s K  TH.-N  T LS  " 

b;  pin 

FOP  J : = 1  UNTIL  Ni  MO 

r  :gin 

3NC ;-N0;  •NDr.-vD;'4;J,:=5(JtJ1'-'MUL  F<K*J); 

PI-PcORT; 

"ND  SWITCH; 

COMMENT  GU  S?  :D  -UMFTI^ns  r/n  u,.5>- ♦ 

FNP  J :  =0  UNTC  M2  DU 
FOR  K  :  =  N  2  (  J  I  UNTIL  M3(J)  DO 
BEGIN 

F<J,KI  :-2.  MXIJN  2+Y<  KU  '2  -  X  (  J  ) -Y  <  K  ))  ; 

DC,  )*K):  =  (X(J)-X(J»  •  2  )  (  Y  (  K  )  -  v  (  K  )  V  V.  ? )  ; 

0  <  2  ,  J ,  K  ) :  =  X  <  J  )  v  ( kf  )  D  ( 1 ,  J  » K  > ; 

G(3,J,K)  :  =  X<  J)  Y <  ^  )  D  (  2  ,  J  , K  )  ; 

D1  l.J.K  ):=(!. -2.  x  <  J  >  )  (Y|K)-Y(n  ■•21; 
ri  R*  X(  JJ-3.  -  X(J)-.  *2)  (Y(K  )  -2-Y(K  C  3); 

1  ;  :  c  ?:-?!■ J^«fr 

rpMM-fT  SOI  V :  =PR  JMT.  G«U  LS  NOW; 
r-PF  I :  =  1  iinjti  M  do 
BEGIN 
DOG  :=C. (  ; 

F?F  J:=l  UNTIL  M?  DP 
=  QF  K  :  =  .M2  (  J  )  +1  UNT  I  L  N 3  <  J  )  DO 
DOG  J  =DOC-  (  =  (  J  ,  K 1  -  0  (  I  ,  J  ,  K  )  -HU  ; 
c  ( T  ) :  =  DC  G ; 

■■'NO; 

e0R  I s - 1  UNTIL  N  DO 
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OI*l< 


FOR  J: =1  UNTIL  N  DC 
BEGIN 

FOR  X:  =  l  UNTIL  M2  DO 

*=OR  L :sN 2 » K  )  +  1  UNTIL  N3( K )  DO 

B(XgJl:*B(XtJ>+(Dl (!«KtL  >*D1( J#K,L )+02( I,K,L  )  'D21  J«fr 

7ND } 

SWITCH; 

ND;«nd; 


i 

i 
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COMPUTER  PROGRAM  5. 
RAYLEIGH-RITZ  METHOD  APPLIED  TO 
V2U  ■  -A2U,  EQ.  (1.1) 


BEGIN 

COMMENT  SOLVES  EIGENVALUE  PROBLEMS  BY  RAYLEIGH-RITZ 
METHOD  USING  TRED2  AND  TGL2  FROM  NUMERI SCHE 
MATHEMATI K  11,  PP.  181-195  AND  PP.  293-306(1968)  BY  MARTIN 
RE INSCH»BOWOLER .HILARY,  AND  WILKINSON. 

FORMAL  PARAMETERS: 


USED 


NUMBER  OF  EIGENVALUES 
AND  THE  NUMBER  OF  EQUATIONS 
MESH  OF  X  VARIABLE 
MESH  OF  Y  VARIABLE 
NUMBER  OF  X  POSITIONS 
NUMBER  OF  Y  POSITIONS 
UNKNOWN  EIGENFUNCTIONS 
ESTIMATING  FUNCTIONS  USED 
COEF.  OF  ESTIMATING  FUNCTIONS 
EIGENVALUES  OF  MATRIX  A 
EIGENVALUES  OF  MATRIX  B 
EIGENVECTORS  OF  MATRIX  B 
EIGENVECTORS  OF  MATRIX  E=T'*A*T 
USED  I Nr TRANS FORMAT  ION  C=T'*D; 

NX  »NY,N»  NPX, NPY! 


FORM  USING 


N 

HX 
HY 
NPX 
NPY 
U(  I  ) 

C  ( I ) 

R  ( I  ) 

02  (  I) 

22(1) 

D<  I  ) 

INTEGER 

REAL  HX,HY , EPS, DOG; 

LOGICAL  T; 

PROCEDURE  TRED2< INTEGER  VALUE  N?  RE AL  ARRAY  D,E(*) 

REAL  ARRAY  A, L ( *, * ) :REAL  VALUE  PU ) I 
COMMENT  REDUCES  SYMMETRIC  MATRIX  TO  TRIDIG. 
HOUSEHOLDER , S  REDUCTION. 

FORMAL  PARAMETERS: 

N  ORDER  OF  SYMMETRIC  MATRIX  A 
D  DIAGONAL  OF  RESULTS 
E  SUB-D I  AGONAL  OF  RESULTS; 

BEGIN 

INTEGER  I , J,K,L; 

REAL  F , G, H , HH , TOL ; 

TOL : =PU/EP  SI  LON : 

FOR  I :=1  UNTIL  N  DO 
FOR  J :  =  1  UNTIL  I  DO  Z (  I  ,  J  )  :  =A (  I  ,  J  )  5 
FOR  1 1 : =N  STEP  -1  UNTIL  2  DO 
BEGIN  I :  =  1 1 ; 

L:-I-2;F:=Z(I ,1-1) :G:=0.0; 

FOR  K : =1  UNTIL  L  DC  G:=G+Z( I , K)*Z( I ,K) J 
H:s=g+F*F; 

IF  G  <=  TOL  THEN 

BEGIN  E ( I ) : =F  : H: =Q. 0 . GO  TO  SKIP;END: 

L  :=L+ 1 ; 

G :  =E ( I ) :  =  I F  F>=  0.0  THEN  -SQRT(H)  ELSE  SQRT(H; 
H:=H-F*G; Z ( I, 1-1) :sF-G:F:=0.0; 

FOR  J: si  UNTIL  L  DO 
BEGIN 

Z  <  J ,  I ) :=Z ( I , J ) /HtGtsO.O: 

FOR  K : =1  UNTIL  J  DO  G : =G+Z < J , K ) *Z ( I , K ) ; 
FOR  K  :  sj+ 1  UNTIL  L  DO  G: =G+Z ( K, J ) *Z < I ,K ) ; 

E ( J ) :=G/H:F:sp+G*Z( J, 1 ) : 

END: 

HH: =F /( H+H ) : 

FOR  j:=l  UNTIL  L  DO 
BEGIN 

F:=Z(I,J);G:=E(J) :=E< J)~HH*F; 

FOR  K :=1  UNTIL  J  DO 

Z< J,K) :=Z( J ,K)-F*E<K)-G*Z< I,K) ; 

END: 

SKIPiDd  ) :  =H  : 


G:=G+Z(I*K)#Z(K»J): 

Z(K»J)is£(K*J)“G*Z(K»l) 


FOR  I*. si  UNTIL  N  DO 
BEGIN  Lssl-l? 

IF  0(1)  -=  0.0  THEN 
FOR  J:  =  l  UNTIL  L  00 
BEGIN  G:=O.G; 

FOR  K : =1  UNTIL  L  DO 
FOR  K : =1  UNTIL  L  DO 
END? 

om:=Z<I,n;z(l, !>:  =  !. 0; 

FOR  J  *  =  1  UNTIL  L  DO  Z ( I  ♦  J ) : =Z ( J , I ) i =0.0 ; 

END; 

END  TRED2  5 

PROCEDURE  TQL2 ( I NTEGER  VALUE  N;REAL  ARRAY  D#E(*); 

REAL  ARRAY  Z(  =>»*)  {LOGICAL  RESULT  FAIL): 

COMMENT  FINDS  EIGENVALUES  AND  EIGENVECTORS. 

FCRMAL  PARAMETERS: 

N  ORDER  TRIDIAGONAL  MATRIX 
D  01  AGONAL  ELEMENTS 
E  SUB-DIAGONAL  ELEMENTS 
Z  MATRIX  OF  HOUSEHOLDER  TRANSFORMATION: 

BEGIN 

INTEGER  I,J.K,LfM: 

REAL  B,C,F,G,H,P,R,S: 

FAIL:=FALSE: 

FOR  I:  =2  UNTIL  N  DO  E  (  1-1 )  :,  =  E ( I ) ; 

E(N) : =B :=F :=0.0 ;  i  , 

FOR  LI : =1  UNTIL  N  DO  i 

BEGIN  L:=Ll:J:=0: 
H:=EPSIL0N*(AB${D(L5)+ABS(E<L) ) ): 

IF  8<H  THEN  8  :  =H ; 

COMMENT  LOOK  FOR  SMALL  SUB-DIAGONAL  ELEMENT; 
FOR  Ml  :  =L  UNTIL  N  DO 
BEGIN  M:=M1 * 

IF  ABS( E( M) )<=B  THEN  GO  TO  CONTJEND; 
CONTUF  M=L  THEN  GO  TO  ROOT; 

NEXTIT : 1 F  J=30  THEN 

BEGIN  FA1L:=TRUE;G0  TO  FIN; END: 

J:=J+1 ; 

P:=(0( L+l)-D( L) )/ (2.*E(L) ) : 

R*-SQRT(P**2+1.0) ; 

H:=D(L)-E(L)/(IF-  P<0.0  THEN  P-R  ELSE  P+R); 
FOR  I :  =L  UNTIL  N  DO  D ( l ) : =Dt I ) -H; 

p  «  « 

COMMENT  QL  TRANSFORMATION; 

P :  =0 ( M) ; C :  =  1. 0 : S; =0.C  ; 

FOR  I:=M-1  STEP  -1  UNTIL  L  DO 
BEGIN 

G J-C*E( I ) :h:=c*p; 

IF  ABS(P)  >=  ABS (Ell) )  THEN 
BEGIN 

C s  =E( I ) /P  :R:  =  SQRT (C#* 2+1.0) ; 

E(  1+1 )  :=S*p*R:$:=C/R:C?=1.0/R: 

END  ELSE 
BEGIN 

C:sp/E( I )  ;R:  =  SQRT (C**2+1.0) ; 

E( l+l )  :  =  S*E(I )*RSS:-1.0/R:Cs=C/RS 
END: 

p;=c*d(I)-$*g: 

D ( I  +  1 1 :=H+S*(C*G+S*D( I)  ); 

COMMENT  FORM  VECTOR; 

FOR  K:=l  UNTIL  N  DO 
BEGIN 

H :  =  Z(  Kf  I+-1)  :Z(K,I  +  1);=S*Z(K«I )  +C*  H: 

2 ( K  ♦  I ) : =C*Z(K i I )~S*H: 

end;end; 

E(L)  :=S*P;D(L  )  :=C*P: 

IF  ABS ( E ( L ) )  >  6  THEN  GO  TO  NEXTIT: 

ROOT :0(L) : =  0 ( L ) +  F ; 

ENO: 

COMMENT  ORDER  EIGENVALUES  AND  EIGENVECTORS: 


FOR  I  *«1  UNTIL  N  DO_t 
BEGIN  K)«nP5«D(I)  ; 

FOR  UNTIL  N  DO 

IF  D(JJ  <  P  THEN  ,  . 

BEGIN  K:*J:P:  «D(J) »END: 
IF  K  -»»  I  THEN 


BtK I  <  I ):D(I):*P; 

FOR  J « *1  UNTIL  N  DO 

BEGIN  P:-Z(JiI) «Z(J»I):=Z(J»K):Z(J»K):=P;END: 
END: END* 

FIN: 

END  TQL2? 

START:READ(N»NPX»NPY»NXtNYtHX*HY  »EPS) : 

BEGIN 

COMMENT  PROGRAM  STARTS  HERE: 

REAL  ARRAY  D1 *D2»E1,E2(1: *NI S 
REAL  ARRAY  Dt R(I : :N» I : :N) : 

REAL  ARRAY  U(1::N»0: : NPX  §0 : : NPY ) ; 

REAL  ARRAY  Et A,B»P( 1 : :N, 1 : :N) : 

REAL  ARRAY  xfo::NPX) : 

REAL  ARRAY  Y(0: :NPY) ; 

REAL  ARRAY  C . Cl. C2 ( I : : N* 0 : : NPX ,0 : : NPY ) : 

REAL  ARRAY  Z1 , Z2 ( 1 : : N » 1 : : N ) ; 

PROCEDURE  FCN: 

COMMENT  COMPUTES  FUNCTIONS  AND  MATRICES  A  AND  B; 

BEGIN 

X(0):=Y(0) :*0.0: 

FOR  I :  *1  UNTIL  NPX  DO  X(I):-I*HX: 

FOR  I :*1  UNTIL  NPY  DO  Y(I):=I*HY: 

FOR  I : =1  UNTIL  N  DO 

FOR  J : =1  UNTIL  N  DO  A ( I , J ) :*B( I , J ) : =0.0: 

FOR  I :=0  UNTIL  NPX  DO 
FOR  J:=0  UNTIL  NPY  DO 
BEGIN 

C(ItI»J):*(X(I)-X(I )**2>*<Y{ J)-Y( J)**2) : 

C(2»I»J):  =  (0. 5-X(  IM*Cfl.  I*  J) : 

CI3.I  ,J)  :  =  <0.f>-Y<  J)  )*C!1,I,JI5 

:  =  d.-2.*X(I  )  )*<Y<  J)-Y(  J)**2); 

C1(2.I ,J)  :  =  (Y<  J)-Y(J  )**2)*<0.5-3.*xm+3.*X(I)**2) ; 

Ci  3,I,J  :s(i.-2.*xiin*(a,5-Y(j)  )*iyijk(ji**2): 

C2( It  I ,J) :=<1.-2.*Y< J )  )*<  X( I )-X( I )**2) : 

C2(2,I,J) Jaj0.5-X(  I)  )*(X(  n-X(n**2»*U.-2.*Y(J) ) : 

C2( 3* I « J ) :  =  (X(I )-X(I)**2)*(0.5-3**Y( J)+3.*Y( J  )**2) : 

END « 

FOR  K:=l  UNTIL  N  DO 
FOR  L : =1  UNTIL  N  00 
BEGIN 

FOR  I : =1  UNTIL  NPX  DO 
FOR  J:=l  UNTIL  NPY  DO 
BEGIN 

A(K*L)  :=A(K.L)+CHK,I  ,J)*C1(L,I,  J)*HX*HY 
+C2 ( K  » I » J  j  *C2 I L  » I,J)*HX*HY: 

BCKjU :=B(K,L)+C(K,I, J)*CU» I, J)*HX*HY: 

END  ♦ 

„  end: 

END  FCN: 

PROCEDURE  WRITER: 

COMMENT  WRITES  ALL  DATA: 

BEGIN 

FOR  I : *1  UNTIL  N  DO 
FOR  J:*l  UNTIL  I  DO 

for8EGIN  DOG:-Zl(I^JI:Zl(I,J):«Zl{J,n:Zl(J,n:*COG:END: 

FOR  I i *1  UNTIL  N  00  D <  L » I > :*Z2( I *L) ; 

FOR  LJ-1  UNTIL  N  DO 
FOR  I  * - 1  UNTIL  N  DO 
BEGIN  DOG: ■0« 0: 

£9?  UNflL  N  00  DOG:*DOG+  ZI(I»J)*D(L»J): 

RU  * L )  :*DOG: 

END: 

FOR  I :■!  UNTIL  N  00 


OMEGA* " «  SQRT(D2(K) ) ) 


BEGIN 

FOR  j:=0  UNTIL  NPX  DO 
FOR  UNTIL  NPV  DO 

BEGIN  U  ( IfJ*L) :  *0 . 0 ; 

wTIL  N  00  U(I»J'LI  :=U(I,J,LI+MI,K)* 
6nd:eno!K,j,u' 

FOR  K:=I  UNTIL  N  DO 
BEGIN  I OCONTROL ( 3)  ; 
tjoIrPr.!  *'l  ’WRITE!  «  ") {WRITE!"  "J; 

w|  T1  .  ")  ;  WR I  TE(  *■  .,?air¥I?9  S^JPUT  5.  (  RAYIE  . . 

KI  .,.«,tII8TOI^:Ms* 

WRITE!"  THE  VALUES  OF  0 ft  I  «  I  - 

WRITE!"  ") JWRITE!"  «);  °F  C(I) 

FOR  Lssl  UNTIL  N  DO 

BEGIN(mTflELDSIZ£:,iL 

”,;"R'rE'"  ") :end:wiute<”  ■•>: 

WRITE!"  « I  s  WRIT*=C  "  »)  ;  U!X,YP'»; 

F°g  ,J:f£  NX  UNTIL  NPX  DO 

BEGIN  INTFIELDSIZE!=2l  D° 

2 NPY.....r*ifeli,ilpOili:;:i!.?f" NM'" 

ENDiENO: 

60  TO  Q; 

END  WRITER;  ; 

FCN ; 

IF  BEGINN  BEGIN  WRITE! "FAIL=1") ; GO  TO  Q; END  ELSE 

™  11:1  BSSTit  fj  gg  01(1!  !=SQRTID1UI  I! 

FOR?J;=l  UNTIL  N  DO  Z 1 ( J , I ) : *21! J , I) /D1 < I ) ; 

FOR  l:«l  UNTIL  N  DO 
FOR  K:=l  UNTIL  N  00 
BEGIN  P( I  ,K) ; =  0.0 J 

END;J:=1  UNTIL  N  00  F<  I  t  K )  :  =  P(I  ,K)+A(  I ,  Jl^l!  J,K) ; 

FOR  I  :  =1  UNTIL  N  DO 
FOR  JJ=1  UNTIL  I  DO 

for  i?li  untilZn(d6J)  :Z1(I,J):=21(J,n;ZltJtns  =  DCGJEND! 

FOR  K J  =1  UNTIL  N  00 
BEGIN  E(I,K):=0,0; 

END  .  J :  =  1  UNTi>-  N  DO  E  !  I  »  K ) 

WRITER*6*  8E<"IN  WRITEt"FAIL-2")  JGO  TO  Q ;  END ; 

END;  * 

Q:GG  TO  START; 

END, 


=E( l »K)+Z1( I , J)*P( j,k> ; 
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COMPUTER  PROGRAM  6. 


DYNAMIC  PROGRAMMING  APPLIED  TO 

(1)  V2U  *  2(x2+y2-X-y)  EQ.  (1.5) 

(2)  V4U  =  X2U  EQ.  (1.3) 


BEGIN 

COMMENT  SOLVES  PARTIAL  DIFFERENTIAL 
EIGENVALUE  PROBLEMS  ON  A  RECTANGLE 
PROGRAMMING  AND  STODOLA'S  METHOD. 
FORMAL  PARAMETERS : 


EQUATIONS  AND 
BY  DYNAMIC 


CM 

OM2 

N 

M 

K1 

VI 


T1 


OMEGA 
OMEGA  SQUARED 
NUMBER  OF  X  POSITIONS 
NUMBER  OF  Y  POSITIONS 
COEFF  IF  NEEDED 
ORDER  OF  THE  OPERATOR 
V1=0  FIRST  ORDER  EIGENVALUE  PRB 
V  1=2  SECOND  ORDER  EIGENVALUE  PRB 
TYPE  OF  PROGRAM  RUN 
T1=0  FIRST  ORDER 
T 1-2  SECOND  ORDER 
Tl-50  EIGENVALUE 

SECOND  NORMAL  DERIVATIVE  ON  BDRY 
FUNCTION  TO  SATISFY 
UNKNOWN  FUNCTION 
NUMBER  OF  EIGENVALUES* 

,ps; 

ERA.OM.OM2,OM02,H,K1; 


5000.0: 


pin 
61 
u 

Cl 

REAL  KE 

REAL  E  ,  . 

REAL  CAT1,  CAT2.C1.C2.C3.D0GI.D0G 2*. 
REAL  ZNORM.ZS; 

INTEGER  N.M.N1.S.T1.V1 »M1; 

INTEGER  01.Q2J  . 

LOGICAL  TOGGLa  5 
REAL  ZNORMl; 

X  J  R  HAD (N.M.H.Kl.Tl.Ql.Vl ) ; 

TOGGLE :=TRUEi 
02: =1; 

CAT1  s  =C,6T2  i  = 

Ml s  s2*  M-2 ; 

0M2: =1000.0; 

N1 s  *0  5 
ERA  :s0.03; 

ERA t  =0  .002; 

6:^0.000001? 

BEGIN 

REAL  ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
I :  =0 
l:*C 
I  :  =0 
I:  =0 
I  :  =0 
J :  =0 
(T1>0 ) 


REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
REAL 
FOR 
FOR 
FOR 
FOR 
FOR 
FOR 
IF 


Mi  ; 


Xl(  0: s N ) ; 

Y 1  ( 0 : :  M ) ; 

U5 ( 0: :N,0 : : M) ; 

U 1 .  U2  ( C:  :N ,  0:  :  M ) ; 

U3  ( C :  :  N .  C  : :  M )  ; 

Cl  .61  ,P1(  Os  *-N.O: : 

Z(0::N+4,0:  :M+4)  ; 

HI  CO: :M»0  ::  M» ; 

11, 0(  1 : : M-l » 1 : :M-1) 
U,B,R(0;  :N,C: : M); 

A.  (  1 :  :N,l:  :  M. l  :  :  M)  ; 
D<l::N.O: :m.q: :mi»; 
CAT (0 : : Mi ; 

C  (  0  : :  N .  0 : 


PI. P3  <0 
P2 , P4 <0 


UNTI  L 
UNTIL 
UNTI  L 
UNTIL 
UNTIL 
UNTI  L 
AND 


: M) ; 

M  )  ; 

Nl  ; 

RE  AD  0  M < U(O.II); 
RE  ADON ( U ( 1,011; 
RE  AD ON ( U( N* I ) ) ; 
READONl  !J<  I  .M)  )  ; 
XI  ( I)  s  =  I*H; 

Y1 ( J )  :  =  J*H  ; 

<  Vl  =  2 )  TH£N 


M 

N 

M 

N 

N 

M 


DO 

DO 

DO 

CP 

DO 

DO 
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OMSGA=".OM) ; 


BEGIN 

FOR  a "Q  UNTIL  M  00 
PI! 1 1 :*0.0; 

FOR  ls«0  UNTIL  M  DO 
P3U  s— Pl(I); 

FOR  I:*0  UNTIL  N  00 
P2( I  :*0.0; 

FOR  I  s=0  UNTIL  N  DO 

P4m:*-P2(n; 

end; 

FOR  Js  *0  UNTIL  M  DO  U1 (I , J ) : =U2( I« J ) : =U3! I , J ) : =U5< I , J ) : *0. 0; 
S:=N; 

BEGIN 

PROCEDURE  WRITER; 

BEGIN  WRITE( "SINGULAR") ; 

GO  TO  W; 

END  WRITER; 

PROCEDURE  WRITl; 

BEGIN 

IOCCNTPOLO); 

WRITE!"  "); WRITE!"  "MWRITE!"  ")  ; 

WRITE!"  COMPUTER  OUTPUT  9. ( DYNAMIC )") ; 

WRITS!"  "); WRITE!"  ");WRITE!"  "); 

INTFISLDSIZE:=l; 

WRITE!"  EIGENFUNCTION*", 02-1 ) ;WRIT£!"  "); 

I NTFIE  LDSI  ZE  :  =2  ; 

WRITE!"  PATH*" »N1,"  OMSGA=".OM); 

WRITE!"  ");WRITE!"  «); 

WRITE!"  X  Y  ", 

"  U !  X  *  Y I " ) ; 

WRITE!"  "); WRITE!"  « ) ; 

FOR  I : =1  STEP  2  UNTIL  N-l  DO 
FOR  J: =1  STEP  2  UNTIL  M-l  DO 
BEGIN  INTFIELDSIZS:=2; 

WRITE!"  ",I  DIV  N * " • " , I  REM  N," 

J  DIV  M,".", J  REM  M ,U  !  I , J )  )  ; 

WRITE!"  ")  ; 

END; 

IF  Q2  =  01+l  THEN  >0  TO  W  ELSE  STARP  ; 

END  WPITl; 

PROCEDURE  STARP; 

BEGIN 

IF  Tl=50  THEN 
BEGIN 

IF  02=1  THEN 

BEGIN  Q2: =2; 

P0R  I :=1  UNTIL  N-l  DO 
FOR  J : =1  UNTIL  M-l  DO 

U!I  ,J)  :=X1!I)*Y1!  J)*!  1.-X1U))*  (1.-Y1  (J)); 

GO  TO  Z5; 

END; 

IF  02=2  THEN 
BEGIN 

FOR  I : =1  UNTIL  N-l  DO 
FOR  J:=l  UNTIL  M-l  DO 
BEGIN 

U2! I , J) :=U( I , J) ; 

U5! I, J ) :=0.0j 

U(I,J):=X1  !ir(l.-Xl(II)*(.5-Xl(I))*Yl!J)«!  i.-Yl(J))  ; 
END; 

CATl: =K f • Q2 := 3 ;Nl:=l;0M2: =5000.0; 

GO  TO  Z5; 

END; 

IF  02=3  THEN 
BEGIN 

FOR  I :  =1  UNT!  |.  N-l  DO 
FOR  j: =1  UNTIL  M-l  DO 
BEGIN 

U  3 !  If J) s  *U( I , J)  ; 

U5( I , J i s  =0.0; 

U!I  ,J)  :*Y1(  JJM1.-Y1  (  J)  )*(  .5-Y1  (  J))*X1  ( I)  t(  1.-X1 !  I) ) ; 
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CAT2?*K^  •02: *4; 0M2:* 5000.0; 

GO  TO  25; 

SND  ;SND ; 

25: 

IF  02>5  THFKJ  GO  TO  Z6;  ... 

COMMENT  PUT  FUNCTION  *0  SAT.  HERE; 

IF  TI-.»50  THEN 
BEGIN 

FOR  I:*0  UNTIL  N  00 
FOR  J: *0_UNTIL  M  DO 

e!?I?J):  =  2.*<  XI  (I  )**2+Yl(J)**2-Xl<n-Yl<Jn; 
C(ItJ):=2.f'*H*H*El(!»  J)  ; 

5ND;END; 

FOR  I s  =  1  UNTIL  M-l  DO 
FOR  J: =1  UNTIL  M-l  00 

M 

IF  }=J  THEN  I1(I.J):-1.0  ELSE  1 1 ( I  •  J ) : =0. 0? 

Q(I  ,J):=C.O; 

IF  f=J  THEN  Q<I » J ) :=K 1*3.0{ 

IF  ABS( I-J  )  =  1  THEN  0(  I « J  )  :=-Kl ; 

END; 

FOP.  J:  =1  UNTIL  N  DO 
FOR  I: =1  UNTIL  M-l  DO 
BEGIN 

IF^i=i*THFN *R ( J  « I  ) : =- 2 .0*K1*U{ J  »0 ) » 

IF  I  =  M-1  THEN  R< J ,1  ): =-2.0*Kl*U< J»M ) ; 

END; 

FOR  I : =1  UNTIL  M-l  DO 
FOR  J: =1  UNTIL  M-l  DO 
A(Nt  I » J) :  =  I 1( I tJ) ? 

FOR  I : =1  UNTIL  M-l  DO 
8 ( N « 1 )  :=-2.0*U(N,I ); 

SWT  TCH * 

Z6 : IF  62  =2  THEN  SWITCH  ELSE  OOP; 

END  STARP; 

PROCEDURE  OOP; 

COMMENT  NORMALIZES  U(X,Y); 

BEGIN 

ZNORM: =0 .0 ; 

FOR  I : =0  UNTIL  N  DO 
FOP  J : =0  UNTIL  M  DO 
BEGIN 

ZS:=ABS(U( I» J ) > ; 

IF  ZS>ZNCRM  THEN  ZNORM:=ZS; 

END; 

FOR  I : =0  UNTIL  N  DO 

FOP  J:  =0  UNTIL  M  DO  U<  I  ,  J)  :  *IJ<  I  ,  JJ  /  ZNORM; 

0MC2 :-QM2 ; 

0M2 :=1 .O/ZMOPM; 

OM: «S0RT(0M2) ; 

IF  Nl>30  THEN 
BEGIN 

WP I  T.E  (  "TOO  MANY  STEPS”); 

GO  TO  w; 

END; 

ZNORMl:=O.Oj 

FOP  I :  =1  UNTIL  N-l  DO 

FOR  J: *1  UNTIL  M-l  DO 

BEGIN  .  ,  . 


ZS:-ABS(Un,J)-U5(ItJ)); 
IF  ZS>ZN0RM1  THEN  ZN0RM1 ; 
END; 

ZN0RM1<ERA  th?N 
BEGIN  KFj-O.O; 

POR  15=1  UNTIL  N  00 

FOR  JJ«1  UNTIL  m  DO  KS:*K 


FOR  J 5  =1  UNTi 

T4]u 

NlJ-Ni+l; 


DO  K5 : *KE+ (U ( I » J )*H)**2; 
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FOP  I :  -0  UNTIL  N  DO 
FOR  J: =0  UNTIL  M  00 
BEGIN 

U5<I ,J) *aU( 1  *  J ) ; 

IF  V1=0  THEN  51(1, J  )s*-U< I,J) 5 
IF  VI =2  THEN  SI (I ,J):=U(I ,J)i 
C(I.J):  =  2.*H*H*(UU.J); 

.  end; 
figure; 

END  OOP; 

PROCEDURE  FIGURE; 

BEGIN 

COMMENT  SOLVES  FOR  U  AMD  B; 

IF  ( T 1 =0 )  OR  ( CT1-50  )  AMD  (V1=0)J  THEN 
BEGIN 

FOR  L S=N  STEP  -1  UNTIL  2  DO 
FOR  I s=l  UNTIL  M-l  DO 
BEGIN  C AT( I ) : =0. 0  5 
FOR  J  J  =  1  UNTIL  M-l  DO 

CAT  <  I  )  !  =C  AT  ( I  )  +  D(L,I,J)N8IL,J)+R(L-ttJ)+C(L-liJ)l! 

B ( L  —  1  ♦  I )  :=CAT ( I ) ; 

END; 

FOR  I:«l  UNTIL  N-l  00 
BEGIN 

FOR  J: =1  UNTIL  M-l  DO 
BEGIN  C  AT ( J ) : =0.0; 

FOR  L :  =  1  UNTIL  M-l  00 

CAT(  J):=CAT ( J)+D< 1  +  1, J,L  ) 

*(U(I-1  ,L)-(B(I  +  l,L)  +  R(I,L)+C(ItL)  )/2.0); 

U ( I , J  )  i =CAT ( J  )  ; 

END; END; END  ELSE  . 

BEGIN 

IF  TOGGLE  THEN 
BEGIN 

U1(0,0»  *®-Pl(  0)-P2(0)  ; 

U1(N,0):=P3(0)-P2(N); 

U 1  (  0 ,  M )  s  —PI  <  MJ+PMO)  ; 

Ul(  N,  M ) :  =P4(M)+P  3( M ) ; 

FOR  I : =1  UNTIL  M-l  DO 
BEGIN 

U1(0.  I )  :=-cl(  I  )+(  U<  0,I  +  1)-2.T'J(0,I  )  +  U(0, 1-1 }  ; 

U1  (N,I):  =  P3  ( I  >  +  (U(N,I  +  l)-2.-*U(N,  I)+U(N,  1-1)  >/(H*H> ; 
END; 

PQR  I J  =  1  UNTIL  M-l  DO 
RSGIN 

Ul( I ,0) : =-P2 ( I )*(U(I+l«0J-2»*U(I ,0)+U( I-l .0 ) >/<H*H)  ; 
Ui(I,M):  =  P4(I )+(U(I  +  l,M  l-2.*U( I ,M ) +U( I- 1 , M ) )/(H*H) ; 
END; 

FOP  I s  =  1  UNTIL  M-l  DO  B1(N,I) :=-2.*Ul(N.I); 

FOR  J:  =]  UNTIL  N  DO 
FOR  I :=i  UNTIL  M-l  DO 
BEGIN  Sl< J. I ) s=n, 0; 

IP  1  =  1  THEN  R1 (J,  I  ):--2.*Ull J,0 ) ; 

IF  IsM-l  THEN  Rl( J,I) :=-2.vU1(J,M); 

end; 

TOGGLE i -FALSE ; 

END; 

POP.  L:=N  STEP  -1  UNTIL  2  DO 
FOR  l:=l  UNTIL  M-l  CO 
BEGIN  CAT ( I ) : =0.0; 

FOP  J : = 1  UNTIL  M-l  DO 

CAT ( ! )  5  CAT  (I  )  +  D<  L,  I, J)*(B1 (L , J)+R1(L-1,  J 1+C(L-1,  J  )  ) 
Blil-l.  I)  :  =  CAT(  I  )  ; 

END; 

FOR  I  :  =1  UNTIL  N-l  CO 
BEGIN 

FOR  J r -1  UNTIL  M-l  DO 
BEGIN  C A T < J) : -0. 0 ; 

FOR  L  :  =  1  UNTIL  M- l  DO 
CAT  ( J  ) :  =C  AT  ( J  )  +'T  (  1  + 1 »  J  »  L  ) 

*(  UK  I-lfL)  -  (  BKI  +  l  ,  L )  +  R1  (I  » L  )  +C  ( I  «L)  )/2.  ); 

Ul( I, J ): =CAT( J  )  ; 
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5NDSEND5  _ 

^FiR^Ss-^uhTIL0®  DO  CI(  It  J  ) :  =  2«*H*H*yi<  I » J) ; 

FOR  LS-N  STEP  -1  UNTIL  2  00 
FOR  I s* 1  UNTIL  M-l  DO 

BEGIN  CAT  (1):  =0.0;  , 

CAT  (II  s  *CAT  ^  !  T  i-D  Vl  II » J  )  +  (B(Lt  J  )+R(L-lt  J)+CI(L-1,  Jl  >  { 

BIL-1 ,1 ):=CAT(I ); 

FOR  ??Si  UNTIL  N-l  00  1 

F§R!Jsal  UNTIL  M-l  DO 
BEGIN  CAT(J):=0.0? 

FOR  L : =1  UNTIL  M-l  DC 
CAT(J  )  :=CAT(J  )+D(  I+l* J«L)  .  t  t  ,,  t 
*(Un-ltL)-(B(I+ltU+R(I,U+CI(ItL)  1/2.);  . 

U(I  ,J> s=CAT(J);  , 

END; END; END; 

IF  02>2  THEN 

BEGIN  ^ 

Cls=DOGl:=DOG2:=0.0;  , 

FOR  l:=l  UNTIL  N-l  DO 

FOR  J:=l  UNTIL  M-l  DO  Cl :=C1+U( I . J)*H*H;  , 

FOR  I : =1  UNTIL  N-l  DO  .  • 

FOR  J:=l  UNTIL  M-l  DO  U(I t J ) : =U ( 1 1 J )-Cl ; 

FOR  I :  =  1  UNTIL  N-l  ; DO 
FOR  J: =1  UNTIL  M-l  DO 

BEGIN  ,  i 

DOGl:=DOGl+U< I* J)*U2( I *J)*H*H? 

DOG2:=DOG2+UU»  J)  *U3!TiJ)*H*H; 

END; 

C2i =D0G1/CAT1 ; C3S -D0G2/CAT2  * 

FOP.  i:«i  UNTIL  N-l  DO 
FOR  J : =1  UNTIL  M-l  DO 
UU*J»s*U(  I,J  J-C2*U2(ItJI-C3*U3(ItJ»; 

END;  i 

IF  Tl=50  THEN  OOP; 

I0CCNTR0U3); 

WRITE!"  ")  ?WRI  TE  (••  WRITS!" 

WRITS!"  COMPUTER  OUTPUT  6.  <  DYNAMIC  )  "  ) 

WRIT”!"  WRITE!"  ");WRITS("  "» ; 

WRITE! "  X  Y  »t 

••  U( XtY  )" ) ; 

WRI  TEC*  ");  WRITS!" 

FOP  I s  *1  STEP  2  UNTIL  N- 1  DO  ! 

C0R  J : =1  STEP  2  UNTIL  M-l  DO 
BEGIN  I MTFI SLDSI ZE : =2  » 

WRITS!"  ".I  DIV  N » "• " » I  REM  N»"  "* 

J  DIV  M «  "  •  "  ♦  J  REM  M,U!I»J)); 

WRITS!"  "); 

END; 

GO  TO  w;  ■ 

END  FIGURE; 

PROCEDURE  SWITCH; 

BEGIN 

INTEGER  L.K; 

REAL  DOM,AMUL,T; 

FOR  L:=M  STEP  -1  UNTIL  2  DO 
BEGIN 
S  *  —  S'— 1  * 

FOR  :I:=1  UNTIL  M-l  DO 
FOP  J:*l  UNTIL  M-l  DO 

D!LfItJ);=0(ItJ )+A(L,!t J) ; 

FOR  I*.»l  UNTIL  M-l  DO 
FOR  J:  =  l  UNTIL  M- 1  DO 
D  (L » 1 ♦ J+M-l ):*Il!ItJ); 

C0R  K :*1  UNTIL  M-l  DO 
PRGIN 

IF  K-M-l  THEN  r-n  TO  Pj 
FOR  U-K+l  UNTIL  M-l  DO 
BEGIN 
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I  F  ABS  < D(  l » K,  K  )  KABS  ( D(  L »  I »K  )  )  THEN 
BEGIN 

FOR  J:,=  l  UNTIL  2*M-2  DO 
BEGIN 

,  T:«D(L«ItJ)SD(LfltJ)i»0(LtKlJ); 

,  P ( L  *K « J ) : =T; 

eMD;5ND;EMD; 

P  : IF  ABS(D(L,K,K)  )<E  THEN'  WRI  TE.P.  ELSE 
BEGIN 

,  DOM:=D(L,K,K)  ; 

FOR  J  i  =  1  UNTIL  2*M-2  DO  0  < L , K .  J ) :  =D ( L . K, 
,  FOR  I: =1  UNTIL  m-i  DO 

BEGIN 

AMULs  =D (L 1 1  ,K) ; 

IF  I=K  THEN  ELSE 
BEGIN 

FOR  J:  =1  UNTIL  2*M~2  DO 
BEGIN 

0(L»I»J):=D(L, It J)-AMUL*D<L, 
END SE NO  SEND  SEND  S?  NDS 
POP.  I :  =1  'UNTIL  M-l  DO 
FOR  J:=l  UNTIL  M-l  DO 
D(L. I.J ) :=D(L . I. J+M-lJ  ; 

FOR  I : =1  UNTIL  M-l  GO 
cOR  J s  =1  UNTIL  M-l  TO 

A(LrliI»J)i"Il( I*J)“0(L»I»J>{ 

END; 

IF  Tl=  50  THEN  OOP  ELSE  FIGURE; 

END  SWITCH; 

STARP;  i  :  l 

Inp; 

W  i  GO  TO  X; 

Z :£ND« 


JJ/DOM 


k,  J ) ; 
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COMPUTER  PROGRAM  7 


DYNAMIC  PROGRAMMING  APPLIED  TO 


(1)  V2u  «  -X2U  EQ.  (1.1) 

(2)  *  8.0  EQ.  (1.7) 


BEGIN 
COMMENT 


SOLVES  PARTIAL  01  FFFRENTI AL 
EIGENVALUE  PROBLEMS  ON  A  RECTANGLE 
PROGRAMMING  AND  STODOLA’S  METHOD. 
FORMAL  PARAMETERS: 


EQUATIONS  AND 
BY  DYNAMIC 


CM  OMEGA 

CM2  OMEGA  SQUARED 

N  NUMBER  OF  X  POSITIONS 

M  NUMBER  OF  Y  POSITIONS 

K1  COEPF  IF  NEEDED 

VI  ORDER  OF  THE  OPERATOR 

V1=0  FIRST  0RD5R  EIGENVALUE  PRB 
Vl  =  2  SECOND  OR  OcR  EIGENVALUE  PRB 
T1  TYPE  OF  PROGRAM  RUN 

T1=0  FIRST  ORDER 
Tl=2  SECOND  ORDER 
Tl=  50  EIGENVALUE 

PU)  SECOND  NORMAL  DERIVATIVE  ON  BDRY 
El  FUNCTION  TO  SATISFY 

U  UNKNOWN  FUNCTION 

Q1  NUMBER  OF  EIGENVALUES; 

REAL  KE.PCi; 

RcAL  E  «cRA,QM.0M2,0M02,H,Kl ; 

REAL  C ATI ,  CAT2 ,C1 ,C2 ,C3, D0G1 ,D0G2 ; 

REAL  ZNORM.ZS; 

INTEGER  N,M,N1,S»T1,V1,M1; 

INTEGER  01.Q2; 

LOGICAL  TOGGLE; 

REAL  ZNORMl; 

X:R”AD (N»M,H,K1,T1 ,01  ,  VI  ) ; 

TOGGLS:=TRUE; 

02 ; =1 ; 

CAT1:=CAT?:=5000.0; 

Ml;  =  2^  m-2; 

0M2: =1000.0; 

Nl ;  =0; 

ERA:=0.03; 

ERA  :=C .002 ; 

E:=0. 000001; 

BEGIN 

REAL  ARRAY 
ARP  AY 
ARRAY 

ARRAY  ui , U2 (C: :M,C: ; M ) 

ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 

ARRAY  0(1  ::N,C: :M,0 
ARRAY  C A T ( 0 : : M ) ; 

ARRAY  C (0 : :M,0: :M ) 

AR°AY  P1«P3(0::M) 

ARRAY 
=0 
«0 


REAL 

REAL 

REAL 

REAL 

REAL 

RCA  L 

PEAL 

REAL 

REAL 

RRAL 

RCAL 

RcAL 

RfAL 

Rr  A  L 

P.FAL 

FfiR 

FOP 

F0° 

=00 

FOR 

FOR 


XI  ( 0 : :M I ; 

Y 1  (  0 : :  m  ) ; 

U5  (0:  s^’.C;  s  M) 
Ul  ,’J2  (C: :  N , C : 
U3(  0:  :  N ,0 : :  v) 
CI,B1,R1(0:  :N 
Z( 0: ;N+4, o: 

H 1 ( 0 i :M.Os : 
11,0(1: :M-1 *1 
U » 3  ,  R  (  0 :  :  N ,  C : 


N,  1:  :  m,  i 


0 : :  M  ) 
:v  +  4) ; 

m>  ; 

:  :m-i  ) 
:”)  ; 

:  :M)  : 

:  :  Ml ) ; 


=  0 
"0 
a0 
*0 


IF  ( T 1 >0 1 


P  2 » P  4  (  0 : :  N  )  ; 

UNTIL  M  00  READON(U (0,1 )) ; 
UNTI L  N  DO  PE ADON(U( I ,01); 
UNTIL  M  DO  0  E  A  D  0  N  ( U ( N,I)|; 
UNTIL  n  do  re  ADONOJ  (  I ,  V  )  ) ; 
UNTI  L  N  DO  XI  (I  )  :  =  ITH  ; 
UNTIL  M  CO  Yl(J 
AMO  <  VI =2 )  TW  =  N 
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BEGIN 

FOR  lt-0  UNT  L  H  00 
Pi(iij«2,*  tvITi  )-Yi  tn**2); 
FOR  i*«0  UNTIL  M  DO 
P3(IIJ—  Pin  ; 

FOR  I J-0  UNTIL  N  DO 
P2U]l>2**(  XI  (I  »-  XI  ( I)**2)  i 
FOR  ts«0  UNTIL  N  DO 
PAU)*—  P2(!  ; 


)  J-P2III; 

FOR  I s »0  UNTIL  N  DO 

FOR  J:*0  UNTIL  M  DO  U1 ( I « J) s  =U2(  I » J1 : =U3 (I  *  J I :  =U5 ( I » J ) : =0 .0 ; 

S:«NS 

BEGIN 

PROCEDURE  WRITER*. 

BEGIN  WRITE ( "S INGULAR" )  * 

GO  TO  w; 

END  WRITER? 

PROCEDURE  WRIT1; 

BEGIN 

I0C0NTR0LC3); 

WRITE  (  "  "UWRITE!"  ");writec  " ) ; 

WRI TE ! "  COMPUTER  OUTPUT  7 .( DYNAMIC  PM  J 

WRIT?  !"  M  I  ;WRITE(  M  ");WRITEt"  ")  *, 

INTFIclDSIZE:=l; 

WRI TE ( "  EIGENFUNCTI0N=",Q2-1 ) ;WRITE("  ") ; 

INTFIELDSIZE:=2; 

WRI  TEC  PATH=",N1,"  OMEGA=" ,  OM.) ; 

WRITE!"  " ) ;  WR I  TE  (  "  "); 


WRITE!"  X  Y  " 

"  U!X,Y)  ")  ; 

WRITE!"  "1;WR1TE<»  ") ; 

FOR  I : =1  STEP  2  UNTIL  N-l  CO 
FOR  J:=l  STEP  2  UNTIL  M-l  DO 
BEGIN  INTFI EL DSIZE : =2 ; 

WRITE!"  ",  I  DIV  REM  N,» 

J  DIV  M  ,  "  t  J  REM  M,1J(  I  ,  J)  )  ; 

WRITE!"  "); 

END; 

IF  Q2=Q1+1  THEN  GO  TO  W  ELSE  STARO; 

ENO  WRITl; 

PROCEDURE  STARP; 

BEGIN 

IF  Tl=50  THEM 
BEGIN 

IF  02=1  THEN 

BEGIN  02: =2; 

FOR  I  sal  UNTIL  N-l  DO 
FOR  J: =1  UNTIL  M-l  DO 

U!I  ,J)  :=X1U)*Y1(  JCU.-X1  !I))*!1.-Y1!  J)  )  ; 
GO  TO  Z5; 

END; 

•  IF  02=2  THEN 
BEGIN 

FOR  I : =1  UNTIL  N-l  DO 
FOR  J : =1  UNTIL  ^-1  DO 
BEGIN 

U?  <  I , J ) : =U( I , J  ) ; 


U5!  I  *  J)  :=0.0; 
U! I t J ) :=X1 ( Ii- 


*(  l.-Xl!  I ))*<  ,5-XHI  )  )*Y1(  J) *! l.-Yl!  J)  >; 


cnu » 

CAT1 : =KP ; Q2 : =3 ;N1:=1;0M2: =5000.0; 
GO  TO  Z5; 

END; 

IF  02-3  THEN 
BEGIN 

FOR  I jaj  UNTIL  N-l  DO 
FOR  J i a i  UNTIL  M-l  DO 
BEGIN 

U3U  , J) :-U( I , J) ; 


U3U  t  J)  :-U(  I ,  J) ; 
f  I  *  J  M-o.o; 

i,ji!»yi  ucn.-Yiun*!  .5-ykj)  i*xim*( 
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CAT2I SKE }Q2;  =  4;0M  2i  =  5000* 0? 

GO  TO  255 
# END; END; 

COMMENT  PUf^F UNCTION 6 tO  SAT.  HERE* 

IF  T l-»s 50  THEM 
BEGIN 

FOR  I  :=0  UNTIL  N  DO 
FOR  J: =0  UNTIL  M  00 
BEGIN 

cul  J )  !=if&*H*H*El<  I  *  J) ; 

END; end; 

FOR  is* I  U^TIL  M-l  DO 

FOR  J: =1  UNTIL  M-l  DO 

IF°I=J  THEN  1 1 < I , J >  s=l .0  ELSE  Il(I,J)s»0.0 
0(1  ,j  ):=0.0; 

IF  I=J  THEN  Q( ! » J > : =K1*3.0 ; 

IF  A8S(  I- J  )  -1  THEN  Q(  I  i  J)  J--K1; 

END? 

Fna  J t =1  UNTIL  N  00 

FOR  ls=l  UNTIL  M-l  DO 

BEGIN 

TF  J  i  =  l  ‘tHEN’rI  J.I  ):  =-2.0*Kl*U(  J«0) ; 
iF  I=M-1  THEN  p(j»I  ):  =-2.0*Kl*U(  J  »M)  i 
END; 

FOR  I : =1  UNTIL  M-l  D 0 
FOR  J: =1  UNTIL  M-l  DO 
A (Ntl ♦ J) :*I 1  (I » J) ; 

FOR  I:=l  UNTIL  M-l  DO 
8(N,I):=-2.0'U(Nt!l; 
cut  fru  • 

Z65  Ic  ^12=2  THEN  SWITCH  ELSE  OOP; 

END  STARP; 

PROCEDURE  OOP;  ,  *  •- 

COMMENT  NORMALIZE?  J(X,Y); 

BEGIN 

ZNORM; =0.0  * 

FOR  I : =0  UNTIL  N  DO 
FOR  J:  =0  UNTIL  M  DO 
BEGIN 

2S:*ASS (U( I«J ) I? 

IF  ZS>ZNQRM  THEN  ZNORM:=ZS» 

END ; 

FOR  I : =0  UNTIL  N  DO 

FOR  J: =0  UNTIL  M  DO  U U t J ) : =U( I ♦ J ) /ZNORM; 
0M02:  =0M2; 

0M2 :=1 .0/ ZN3RM; 

OM: =$Q  &  T( 0M2 )  ; 

IF  Nl>30  TH=M 
P5GIN 

WRI  T€  (  "TOO  MANY  STEPS'*); 

GO  TO  W; 

CKJp  * 

ZNORMi ; =0.  0; 

FOR  Is  *1  UNTIL  N-l  DO 
FOR  J: =1  UNTIL  M-l  DO 
BEGIN 

ZS:=ABS(U( I.J )-U5 ( I » J ) ) ; 

IF  ZS> ZNORMI  THEN  ZN0RM1:=ZS; 

END ; 

IF  ZN0RM1<ERA  THEN 
BEGIN  KE : =0.0; 

PH Q  t  •  si  UNTIL  N  no 

FOR  J:=l  UNTIL  M  DO  KE: =K c+ ( U ( I t J ) *H)*»2 ; 

WRI Tls 

END; 

N1 : =Nl+l; 
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FOR  I  S -0  UNTIL  N  DO 
FOR  J: =0  UNTIL  M  DO 
BEGIN 

U5<  I  *  J)  2  “U  <  I  *  J ) ; 

IF  Vl^O  THEN  51<I.J):=-U<  I.J)  ; 

IF  VI =2  THEN  Fltl.J >:-U(I,J)i 
C(I  »J)  :=2.*H*H*E1<I  ,J)  ; 

END ; 

FIGURE 5 
END  DOP; 

PROCEDURE  PIGURE; 

BEGIN 

COMMENT  SOLVES  FOR  U  AND  B; 

IF  (Tl-0  )  OR  (  (T 1=50  >  AND  (V1=0))  THEN 
BEGIN 

FOR  Li=N  STEP  -1  UNTIL  2  DO 
FOR  1 :=!  UNTIL  M- 1  DO 
BEGIN  C AT ( I) :=C.C; 

FOR  J:  =  1  UNTI  L  M-l  DO 

CAT(I):=CAT( I )+D<L,I, J>*(B(Lt J  »+R (L-l , J )+C( L-l , J )) ; 
B(  L-l  »I ) i =CAT ( I  l; 

end; 

FOR  1 2=1  UNTIL  N-l  DO 
BEGIN 

FOR  JJ-1  UNTIL  M-l  00 
BFGIN  CAT { J ) : =0 «0 ; 

FOR  LS=1  UNTI L  M-l  DO 

CAT  (Ji:=CAT(JJ+D(l  +  l*J»L) 
*(U(l-l,L)»(S(I+l,L)+R(l,L)+C(nL)l/2.0)i 
U<  T  » J ) s=C  AT ( J 1  ; 

snd;end;2nd  else 

BEGIN 

IF  TOGGLE  THEN 
BEGIN 

UKC.O)  :=-Pl(0)-P2<0)  ; 

Ul(N.0»t-P3(O)-P2(t! ); 

U1(0»M) : =-pi ( m ) +P4 (0 ) ; 

U1 ( N  »  M)  5=P4(N)+°3( 

FOR  I :=1  UNTIL  M-l  DO 
BEGIN 

UUOiHss-PlU)  +  (ij|0|Hl)“2»*UIO»U+UlOtI-lM  /  (Hr*H) ; 
U1(N*I  )  s  =  P3(I)  +  (U(N,  I+1>-2.*U(N,  1  >+U(N,  1-1)1 /  <H*»‘H)  ; 
END ; 

FOR  l:=l  UNTIL  N-l  CO 
BEGIN 

Ul(  I  ,C>  s  =-P2  (t ) +( U< I +1 ,0)-2,*U(I,0  KUU-1,0)  )/<H*H); 
Ul(l.V)s«PAU  »  +  <U(  I*1.N)-2.*UII  .N|+UU-l.Mn/(H+H)  ; 
end; 

fop.  I  :  =1  UNTIL  M-l  DO  Bl<N,I)i  =  -2.*Ul<Ntn; 

FOR  J: = 1  UNTIL  N  DO 
FCR  I : =1  UNTIL  M-l  DO 
BEGIN  Pit  J,I)  :  =0.  0 ; 

IF  1=1  THEN  Rl<J,ll:=-2.*Ul(Jt0>; 

IF  1»M-1  THEM  Rl(  J,H  S=-2.AU1 
END ; 

TOGGLE :  =FAL$E ; 

END; 

FOP  L:=N  STEP  -1  UNTIL  2  DO 
FOR  I J -1  UNTIL  M-l  DO 
BEGIN  C. A T  ( I  )  :  =0.0; 

FOR  J : = 1  UNTIL  M-l  DO 

CAT  (I  )  :=CAT  ( I  )+0(L  t  I.  J  >*(  B1  (L,  J  l+RKL-l,  J  )+C<L-l»  J  )> 
Bl<  L— - 1*1)  :=CAT(!  ) ; 

END; 

FOP  I :=1  UNTIL  N-l  DO 
BEGIN 

FOP  J  s s 1  UNTIL  M-l  DO 
BEGIN  CAT ( J) :=0.0 ; 

FCR  L :  =  1  UNTI  L  M-l  DO 
CAT (J) :=CAT (J )+D(  I  +  l, J»L ) 

t(UlCI-l,L)-(?l(I+l,L)+Rl  (I.U+CHiLI  1/2.); 

U 1 < I  *  J ) :=CAT( J  )  ; 
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6ND;fnd; 

FOB  I:«0  UNTIL  N  00 

FOR  J: =C  UNTIL  M  DO  Cl ( I . J ) :=2.*H*H*U1( I . J) i 
FOR  LtsN  STEP  -I  UNTIL  2  DO 
FOR  I :=1  UNTIL  M-l  DO 
BEGIN  CAT ( I ) t-0 .0  ; 

FOR  J :  =1  UNTIL  N-l  DO 

CAT ( 1 )  i=CAT(  n+D(LTI  •  J  )*<  9( L*  JI  +  R(L-l .  J  >+CI  ( L-3 
B<L-ltn:=CAT(n; 

END; 

FOR  l:«l  UNTIL  N-  1  DO 
BEGIN 

FOR  J:=l  UNTIL  M-l  00 
BEGIN  CAT ( J ) : =0.0; 

FOR  L; =1  UNTIL  M-l  DO 
CAT!  J ) :=C ATI J) +0  < 1  +  1 ♦  J,L> 

♦  (UU-l.L  )-<B(  l+l.U+R!  I,L)+CI  (  IiLII/2.  I 
U<  I  ♦  J ) ;  =C  AT  ( J ) ; 

END  ;END  ;END  ; 

IF  C2>2  THEN 
BEGIN 

C1;=DOG1:=DOG2:=C.O; 

FOR  IJ=1  UNTIL  N-l  DO 

FOR  J ;  =  1  UNTIL  M-l  DO  Cl  :=C1 +U ( I , J >*H*H; 


FOR  l:=l  UNTIL  N-l  DO 
FOR  J :  -1  UNTIL  M-l  DO  U  (1  ,  J  i :  «U  ( I  »'J  )-Cl  5 

“  '  -  "  '  DO 

DO 


FOR  I s-1  UNTIL  N-l 
FOR  J; =1  UNTIL  M-l 
BEGIN 

D0G1 :  sDOG l  +  UI  I  ,  J)*U2U  ♦ 

D0G2 :  =D0G2+U  (  I ,  J  )  *U3 1  I ,  J  J  *H«-H  ; 

END; 

C2; =D0G1/CAT1 ;C3: =DCG  2/CAT2 ♦ 

FOR  1 5*1  UNTIL  N-l  00 
FOR  J:  =  l  UNTIL  M-l  DO 
U(  I  r  J  )  S  -U(  I  *  J  1  -C  2M  U2(  I  ,J)-C3*U3(I  .J); 

END; 

IF  Tl- 50  THEN  OOP; 

ICCCNTROL ( 3 ) ; 

WRITE!"  “  I ;  WP ITS!"  ");VJRITC(«  «); 

WRITE!  *•  COMPUTER  OUTPUT  8 

WRITS!"  " ) ; W R I T g { "  " ) ; WRITE!"  «); 

WRITE!"  X  Y 

"  U( X»Y)") ; 

WRITE!"  ") ;WRITE("  " ) 5 
FOR  I : =1  STEP  2  UNTIL  N-l  DO 
FGR  J:=l  STEP  2  UNTIL  M-l  DO 
BEGIN  INTFIELOSUc:  =2; 

WRITE!"  "tl  DIV  N,".",I  REM 

J  DIV  «,«,«,  J  REM  M,U(  ItJli  ; 

WRI  TE  (  "  ")  ; 

END ; 

GO  TO  W; 

END  FIGUR-; 

PROCEDURE  SWITCH; 

BEGIN 

INTEGER  L » K  i 
REAL  DOM.AMUL.T; 

FDR  L:=N  STEP  -1  UNTIL  2  DO 


BEGIN 

S ;  -  s—  1  * 

FOR  I : =1  UNTIL  M-l  DO 
FOP  J : =1  UNTIL  M-l  DO 

D (L» 1  * J ) :=0 ( I  * J )  +  A(L»I.  J) ; 
FOR  I : =1  UNTIL  M-l  DO 
FOP  Js=l  UNTIL  M-l  DO 
D(L.I.J+M-1):  =  I1U,J); 

FOR  K ; =1  UNTIL  M-l  DO 
BEGIN 

IF  k=m-i  THEN  GO  TO  P; 
FOR  1:=K+1  UNTIL  M-l  DO 
BEGIN 


t  J )  ) ; 


.  (DYNAMIC  I") 


N. " 
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IP  ABS(D(L.K«KJ  )<AB S ( D{  L  •  1 • K )  )  THEN 
BEGIN 

FOR  Ji-1  UNTIL  2*M-2  DO 
BEGIN 

T:=D(LtI*J)»D(LiI«J)S=D(L»KfJ); 
0  (  L  ♦  K ,  J)  :=T; 

END*  END ;END  $ 

P:IF  ABS(0(L.K*Kn<«  THEN  WRITER  ELSE 
BhGIN 

DOM:=0(L,K,K); 

E8S  i\:\  until  £W°  D'L’K’J,!=D,L-K' 

BEGIN 

AMUt:  *DIL.I*K); 

IF  I  =  K  THEN  ELSE 
BEGIN 

FOR  J t  sl  UNTIL  2*M-2  DO 
BEGIN 

r.(A  D(L»ItJ):=D(L#l»J)“AMUL*0(L, 

END jEND; END; END; END; 

FOR  It =1  UNTIL  M-l  DO 
FOR  J:  =  l  UNTIL  M-l  DO 
„  D(L,1,J):=D(L,I,J+M-1); 

FOR  X:*l  UNTIL  M-l  DO 
FOP  Jt  =  l  UNTIL  M-l  DO 

AJL-ltl *  J  » s  =  1 1 ( I* J)-0(L,I,J); 

IF  Tl=  $0  THEN  DOP  ELSE  FIGURE; 

END  SWITCH; 

STA RP; 

END; 

end; 

W;G 0  TO  x; 

ZtsND. 


JI/DOMJ 


K,u ) ; 
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COMPUTER  PROGRAM  8. 
GALERKIN’S  METHOD  APPLIED  TO 
72U  »  2(x2+y2-x-y )  EQ.  (1.5) 


FORMAL  PARAMETERS* 


N 
M2 
M3 
M2 
N3 
F 

cm 

mv 


'^W&f'TT^sn.aNs 

SfeivN^BPo«IS  l 

ARRAY  OF  UPPER  Y 
KNOWN  FUNCTION 

C&PI  lO^ Tls?lS.ATriNG°  EuAcT  IONS ! 


POSITIONS 

POSITIONS 

POSITIONS 


EtDOif  »CAT  »0E  VtAMULiTj  DOM* 


REAL  HjE 
REAL  Hi; 

REAL  H2,H3*. 

INTEGER  M2,M3{ 

INTEGER  N,M,R; 

E :=0. OOOOOl ; p :;0; dev: 

READIN  »M2»  M3«H2»H3)» 
Hli -H2*H3» 


l.o; 


BEGIN 

REAL  ARRAY 


REAL 

REAL 

REAL 

REAL 


ARRAY 

ARRAY 

ARRAY 

ARRAY 


X<0::M2»; 

v  <os  j  M3 1 ;  , , 

8(is:M»lJ : N+l I! 

A«C(l::N>? 

FtU(0  ::M2  *0: :  M3); 

2) 


real  ARRAY  O.Dlll:*N*0«; 
*NT EGER  ARRAY t N2«N3lO: :M 


:  M3) 


FOR 

FOR 

FOR 

FOR 

FOR 

FOR 


I  5  =0 
U=0 
i;  =0 
i ;  =0 
I :  =0 
J:  =0 
BEGIN 
P(ItJ) 


UNTIL 
UMTI  L 
UNTIL 
UNTIL 
UNTIL 
UNTIL 

=C.O 


M2 

M3 

M2 

M2 

M2 

M3 


DO 

DO 

DO 

DO 

DO 

DO 


XtI)t=I*H2; 

YU  ) :  si  <H3; 
PEAD0N(N2(  j  }  )  *. 
RFADGNl  N3  1 1 ) )  ♦ 


N  00, ,  „ 

D( L« I *U  l :=Dl(L.i ♦  J ) !=0«0* 
BEGfN^dOMMENT  SOLVE  HERE? 

BEG IN° St?S InSuLAR"  MGC  TO 
IKBcEDJJr^ORM; 

BEGIN  COMMENT  GIT  U ( I , 

WRI  Tc  (  M  •*)  ;  WRITE  ("  M> ; 

WRITE(M 
WRI  TEC’ 

WRITE(« 

writ? ( "  w ) ; 

FOR  15=1  UNTIL  N  00 

MI  1 V  =  B <  I  ♦  N+l ) ; 

INTFISLDSI ZE:-lt 
WRI TE  l" 

FOR "  J  ?=0  UNTIL  M3-1 
FOR1 1 s  s0  UNTIL  M2-1 


«•)*,  WRITE  ("  “)  5 


o; 


J)  AND  A ( I )  5 
WRI  (  w  »• )  ; 

COMPUTER 


OUTPUT  10. (GALSRKIM) M) 


VALUES  OF  Cll)  ARE”) ;  WRI  TE  ( **  ">J 


C".I."*".A(1)) ;WRITE<*  "> 


DO 

DO 
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BEGIN  0OG:*O*Oj 

FOR  L :  =  1  UNTIL  N  DO  DOG:=DOG-t-A(  U*D(  L  .1  .  J> 
lUtJl  J-D0G5 
ENO;END; 

WRITE! "  "  I  ;WR  ITS  ( '•  "  );WRITE!"  "1; 

WRITE!"  X  V 

"  U!  X  •  Y I "  I  ; 

WRITS!**  "I  ;WPITg!*«  *•  )  i 
FOR  I : =1  STEP  2  UNTIL  M2-1 
FOR  J: =1  STEP  2  UNTIL  M3-1 
BEGIN  INTFIELDSIZ5J  =2*- 
WRI  Tc(" 


J 

**)  { 


M3 


DO 

DO 

",I  DIV  M2  , "  •  "  ?  I  R£M  M2  i 


",  I  DIV  M2  , " . " ,  I 
J  REM  M3,U(I ,  J)  I  ; 


B(I,M1  +  1):  =  C(  II 


THEN 


B  ( K ,  J  ) 


DIV 

WRITE!" 
end; 

GO  TO  0; 

END  PERFORM; 

PROCEDURE  SWITCH; 

BEGIN  COMMENT  SOLVES  FOR  A!  I  I  COEFF; 

INTEGER  Nl,Mi; 

Ml:=N;Nls=N+l ; 

FOR  I i =1  UNTIL  Ml  DO 
FOR  Ki -1  UNTIL  Ml  00 

BEGIN  IF  K«M1  THEN  GO  TO  P; 

FOR  It-K+1  UNTIL  Ml  DO 

BEGIN  IF  APS!  B(  K,  K)  KABS!  B(  I  ,K)  ) 

BEGIN  R:-R+l; 

FOR  J : =1  UNTIL  M  DO 

BEGIN  T:=B(I , J) ;B(I  *J) s=B(K,J) 

end;End; 

P: IF  ABS(B!K,K)J<S  THEN  WRITFR  FL SE 

BEGIN  D£ V : ®DE V*B(K,KJ;DCM:=8!K,K); 

FOR  J: =1  UNTIL  N1  00  B ( X , J ) : =B ! K , J ) /DOM ; 

FOR  I  !=1  UNTIL  Ml  00 

BEGIN  AMUL:=B(I,K> ; 

IF  I=K  THEN  ELSE 
BEGIN 

FOR  J:  =  l  UNTIL  N1  DO 

BEGIN  R!I,J)s=B!I,J)-AMUL*B!KtJ) ; 

Hnd;end;end;£ND;snd; 

PERFORM; 

ENC  switch; 

COMMENT  GUESSED  FUNCTIONS  0! I >  HERE; 

FOP  J : =0  UNTIL  M2  00 
FOR  Ks -N2 ( J  J  UNTIL  N3(J)  DG 
BEGIN 

FI J ,K ) S  =2. X( J)**2+Y(K)**2-X( Jl-Y! K) I ; 

0(1 ,  J,K):  =  {  X!  J  >**2-X(  J  )  3 )  *(  Y (  K  )  —  Y (  K  )**2J  ; 

D(  2  ,  J,  K)  S  =<  X!  J)  -X  (  J  )**2  I*  ( Y  (K  ) -*2-Y  (  K  )  **3  )  ; 

D<3.  J.K  ):**(  X(  J  l**2-X<  J)**3)*(  Y(  K)**2-Y(  K)**3)  ; 

D1  <1,J,K) i  =  (2.-6.*X (J)  >*<Y(K}-Y(K  )**2)  + 

2.«<  X!  J  )**3-X<  J  ! ■*  * 2  )  ; 

01  (2. J.K ) S=(2.-6. *Y!K ) )*{ X! Jl-X! J )**2)  + 
2.*(Y(*)**3-Y<K)*'*2)  ? 

D1(3,J,K)  s  =  (  2.-6.  *X<  J)  J:M Y!K I  **2-Y(K  1**31  + 
<2.-6.*Y(K  J )  ->  ( X  ( J  )**2-X(  !}**3); 

END ; 

COMMENT  SOLVE  FOR  INTEGRALS  F*D  AND  LD*D; 

FOR  1**1  UNTIL  N  DO 
BEGIN  DOG : =C«  C; 

FOR  J : = 1  UNTIL  M2  DO 

FOR  K:-N2<J)+1  UNTIL  N3(J)  DO 

DOG  : sDOG+ (  F  <  J , K ) *D ( l , J,K)*H1 ) ; 

C ( I  )  :  =  DOG; 

FOR  J :-l  UNTIL  N  DO 
BEGIN  CAT:=0. 0; 

FOR  K : = 1  UNTIL  M2  DO 

FOR  L:=N2(K)+1  UNTIL  N3(K)  DO 

CAT :=CAT+ ( 0 1( J , K , L )*D ( I . K . L I f HI  I ; 

8  ( I ,  J  I : sCAT ; 

END; 

IF  I=N  THEN  SWITCH; 


t;end 
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